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1.  SUMMARY 


We  developed  a  unified  source  moment  tensor  and  tomographic  inverse  method 
based  on  finite-difference  simulations  of  wave  propagation  at  local  to  teleseismic 
distances.  The  new  method  accounts  for  complex  wave  propagation  in  three-dimensional 
(3D)  earth,  linearizes  the  inverse  problem  by  iteratively  updating  the  earth  model,  and 
provides  an  accurate  way  to  integrate  body-  and  surface-wave  observations.  These 
developments  include: 

(1)  A  non-staggered-grid  finite-difference  code  in  the  spherical  coordinate  system  for 
simulation  of  regional  and  global  wave  propagation  (Zhang  et  al.,  manuscript  in 
preparation,  hereinafter  Zhang  et  al.,  2010), 

(2)  A  new  and  efficient  implementation  of  the  Perfectly  Matched  Layer  (PML)  that 
provides  a  better  performance  in  simulation  of  waves  incident  at  a  near-grazing 
angle,  low  frequency  waves  and  evanescent  waves  (Zhang  and  Shen,  submitted, 
2009;  Appendix  A), 

(3)  Earthquake  location  and  moment  tensor  inversion  using  hierarchical  Finite- 
Difference  Strain  Green  tensor  (FD-SGT)  databases  (Shen  et  al.,  submitted,  2009; 
Appendix  B),  and 

(4)  Full- wave  tomographic  inversion  of  body  and  surface  waves  for  a  joint  P  and  S 
velocity  model.  This  full-wave  tomography  is  based  on  an  accurate  representation 
of  finite-frequency  Pn/Sn  waves  (Zhang  et  al.,  2007),  a  fuller  utilization  of  an 
arrival  on  all  three  components  of  seismic  records  (Shen  et  al.,  2008a),  and  the 
dependency  of  P  waveforms  on  S  velocity  anomalies  (Zhang  and  Shen,  2008; 
Appendix  C)  and  Rayleigh  waveforms  on  P  velocity  anomalies.  As  an  example 
we  have  successfully  applied  the  full-wave  tomography  method  to  the  southeast 
Tibetan  plateau  (Shen  et  al.,  2008b;  Zhang  et  al.,  2009;  Appendix  D). 

The  finite-difference  methods  developed  to  date  make  it  possible  to  construct 
hierarchical  FD-SGT  databases  efficiently  and  with  great  flexibility.  The  spherical  finite- 
difference  method  is  well  suited  for  regional  and  global  wave  propagation,  while  the 
finite-difference  method  with  boundary-confonning  grids  that  follow  topography 
provides  higher  accuracy  at  local  and  near  regional  scales.  The  self-consistent  FD-SGT 
databases  constructed  from  finite-difference  simulations  of  wave  propagation  in  full- 
wave  tomographic  models  can  be  used  to  detennine  the  moment  tensors  within  minutes 
after  a  seismic  event,  making  it  possible  for  real  time  monitoring  using  3D  models. 

2.  INTRODUCTION 

Accurate  3D  earth  models  are  essential  for  accurate  predictions  of  seismic 
observables  and  source  characterization.  While  recovering  the  true  form  and  magnitude 
of  earth  structure  has  always  been  the  goal  of  tomographic  inversion,  recent  studies  show 
that  methods  based  on  full  wavefields  in  3D  reference  models  have  several  distinct 
advantages  in  relating  seismic  observations  to  highly  heterogeneous  Earth  structures. 
Made  possible  by  advance  in  high-perfonnance  computation,  this  full-wavefield 


1 


approach  accounts  for  complex  wave  propagation  in  3D  heterogeneous  earth  (Zhao  et  al., 
2005;  Tromp  et  al.,  2005;  Liu  and  Tromp,  2006;  Zhang  et  al.,  2007),  enables  fuller 
utilization  of  an  arrival  on  all  three  components  of  seismic  records  (Shen  et  al.,  2008a), 
and  linearizes  the  inverse  problem  by  iteratively  updating  the  3D  reference  model.  An 
important  benefit  of  physically  realistic  and  accurate  modeling  of  full  wavefields  in  3D 
models  is  the  consistency  of  the  system  of  equations  in  inversion.  This  is  particularly 
important  for  the  integration  of  different  types  of  observations  (P,  S,  surface  waves 
including  empirical  Green’s  functions  derived  from  ambient  noise)  in  inversion,  which  is 
an  essential  step  to  obtain  a  coherent  and  self-consistent  model  of  the  crust  and  mantle. 
Seismic  tomography  based  on  path-  or  structure-averaging  (rays)  and  ID  or  2D  reference 
models  does  not  accurately  place  the  3D  sensitivities  of  seismic  observations  in  a 
complex  3D  model.  The  mismatch  in  sensitivities  due  to  approximations  in  theory 
becomes  problematic  when  the  target  resolution  is  at  scales  smaller  than  the  Fresnel  zone. 
For  example,  a  1-Hz  P/Pn  wave  recorded  by  a  receiver  -1000  km  from  the  source  has  a 
Fresnel  zone  width  of  -200  km  in  the  middle  of  the  wave  path,  van  der  Lee  and 
Frederiksen  (2005)  recognize  that  a  more  consistent  and  linear  system  of  equations 
offered  by  the  full-wavefield  approach  would  fit  the  combined  data  sets  better  under 
comparable  regularization  biases. 

In  this  field  two  methods  have  been  developed  so  far  to  carry  out  full-wave 
tomography  iteratively  with  3D  reference  models.  One  is  the  adjoint- wavefield  (AW) 
method,  which  back-propagates  the  data  from  receivers  to  image  structure  (Tromp  et  al., 
2005;  Liu  and  Tromp,  2006).  The  other  is  the  Strain  Green  Tensor  (SGT)  based, 
scattering-integral  (SI)  method,  which  calculates  and  stores  the  sensitivity  kernels  for 
each  data  functional  (Zhao  et  al.,  2005;  Zhang  et  al.,  2007).  Both  methods  are  based  on 
the  same  full  wave  theory.  The  main  differences  are  in  their  computational  approaches. 
In  general,  the  adjoint  method  is  CPU  intensive,  while  the  SGT  database  approach 
requires  a  large  disk  space  and  fast  network.  Chen  et  al.  (2007a)  provided  a  comparison 
of  the  AW  and  SI  methods.  We  note  that  an  advantage  of  the  SI  method  is  that  the  SGT 
databases  can  also  be  used  in  source  moment  tensor  inversions,  making  it  possible  to 
account  for  complex  wave  propagation  in  3D  reference  models  in  real  time  seismic 
monitoring.  To  date,  the  AW  and  SI  methods  use  two  different  approaches  to  calculating 
synthetic  waveforms:  the  spectral-element  method  (SEM)  (Komatitsch  et  al.,  2002)  and 
finite-difference  method  (Olsen,  1994;  Zhang  et  al.,  2008;  Zhang  et  al.,  2010), 
respectively. 
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METHOD,  ASSUMPTIONS,  AND  PROCEDURES 

3.1  Overview  of  Unified  Tomography  and  Source  Moment  Tensor 
Inversion  Based  on  FD-SGT  Databases 

In  the  SGT-based  approach,  the  displacement  can  be  expressed  as  (Zhao  et  ah,  2006), 
M„(Grs)  =  ^G„/(r;rs)M.i. 

=  iKG,„(r,;r)  +  ^G,„(r,;r)]M/,  (1) 

=  Ht„(  r,;r)M„ 

where  r  is  the  location  vector,  r  is  the  source  location,  Mjt  the  source  moment  tensor 
component  (i,j=\  -3),  G(r;r)  and  G(r;r)  the  Green  tensors  for  a  unit  impulsive  force  at 
the  source  andr,  respectively,  and  Hijn(rs;r)  the  strain  Green  tensor  composed  of  the 

spatial  gradient  of  the  Green  tensors.  The  subscript  n  denotes  the  component  of  the 
displacement.  The  source-receiver  reciprocity  and  the  symmetry  of  the  moment  tensor 
have  been  applied  in  Equation  1 .  Notice  the  SGT  can  be  at  any  location  in  the  medium 
and  thus  can  also  be  at  the  receiver  location  rR  .  Equation  1  is  the  basis  upon  which  the 
FD-SGT  moment  tensor  inversions  are  developed  (Zhao  et  ah,  2006;  Shen  et  ah,  2009). 

The  wavefonn  change  at  rR  due  to  a  local  model  perturbation  at  an  arbitrary  position 
r  ,  can  be  expressed  as  (Zhao  et  ah,  2006): 

=  Hijn^rR)^Cijki{r)£ik{Y-rs)  (2) 

where  SCjk,( r)  is  the  perturbation  of  the  fourth-order  elasticity  tensor  and%(r;rs)  is 

spatial  gradient  elements  of  the  displacement  vector,  which  requires  one  forward  wave 
simulation  from  the  source.  The  total  waveform  change  due  to  the  model  perturbation  in 
the  entire  volume  in  the  model  can  be  found  by  volume  integration  of  Equation  2.  The 
above  waveform  perturbation  equation  is  the  essential  ingredient  in  the  SGT-based 
tomographic  method.  The  Frechet  kernels  relating  elastic  and  anelastic  properties  to  any 
seismic  data,  such  as  the  frequency-dependent  traveltime  and  amplitude  anomalies,  can 
be  obtained  from  it  (Zhao  et  ah,  2005;  Zhang  et  ah,  2007;  Zhang  and  Shen,  2008). 

Figure  1  provides  a  schematic  illustration  of  our  unified  and  self-consistent  source 
moment  tensor  and  tomography  inversion  based  on  3D  FD-SGT  databases.  Station  SGT 
databases  are  constructed  from  an  initial  3D  reference  model  by  finite-difference 
simulation  of  the  responses  to  orthogonal  unit  impulsive  point  force  acting  at  the  stations. 
The  SGTs  in  a  small  volume  surrounding  the  source  reference  location  extracted  from  the 
databases  make  it  possible  to  invert  for  source  mement  tensors  and  location  in  a  global 
optimization  scheme.  Forward  wave  propagation  simulation  from  the  source  in  the  3D 
reference  model  provides  synthetic  waveforms  at  the  station  and  within  the  entire  model. 
Travel  time,  amplitude,  and  other  waveform  anomalies  are  measured  from  the  observed 
and  synthetic  waveforms  at  stations.  The  forward  wave  field  (from  the  source)  and  the 
station  SGTs  are  used  to  calculate  the  finite-frequency  sensitivities  of  the  measurements 
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to  perturbations  in  Vp,  Vs  (or  bulk  and  shear  moduli),  density,  and  attenuation.  Together 
with  the  sensitivities  to  source  parameters,  the  measurements  and  structural  sensitivity 
kernels  are  used  to  invert  for  the  earth  structure.  The  tomographic  inversion  results  are 
added  to  the  3D  reference  model  and  validated.  This  process  can  be  repeated  to 
progressively  improve  the  model  resolution. 

3.2.  Finite-Difference  Wave  Simulations 

An  accurate,  efficient,  and  easy-to-use  numerical  method  is  the  most  basic  tool  in 
seismic  studies  using  3D  reference  models.  Researchers  have  developed  several  methods 
to  simulate  wave  propagation  in  3D  models.  These  include  finite-difference,  pseudo 
spectral,  spectral  element,  and  finite-element  methods  (e.g.,  Hung  and  Forsyth,  1998;  Igel 
et  ah,  2002;  Komatitsch  et  ah,  2002).  In  particular,  the  spectral  element  method  (SEM) 
has  been  used  in  the  adjoint- wavefield  tomography  and  source  moment  tensor  inversions 
(e.g.,  Liu  et  ah,  2004;  Tromp  et  ah,  2008).  The  SGT  databases  can  be  constructed  using 
any  numerical  method,  though  we  have  focused  on  the  development  of  the  finite- 
difference  method  for  several  reasons.  The  finite-difference  method  does  not  require  the 
grid  mesh  to  conform  to  interior  boundary  topography  (Moczo  et  al.  2007).  It  is  thus 
relatively  straightforward  to  apply  the  method  to  various  regions.  The  finite-difference 
computation  can  be  easily  parallelized  on  high-performance  distribute -memory 
computers.  It  is  also  one  of  the  most  widely  used  numerical  methods  in  seismological 
studies. 

Previous  finite-difference  simulations  of  wave  propagation  in  spherical  earth  were 
based  on  staggered  grids,  which  needed  interpolation  to  generate  values  not  located  on 
staggered  location,  causing  a  loss  of  precision  (Igel  et  al.,  2002).  To  overcome  this 
problem,  we  implement  a  finite-difference  algorithm  with  non-uniform,  non-staggered 
grids  in  a  spherical  coordinate  system  (Zhang  et  al.,  2010).  No  interpolation  is  needed 
during  the  calculation  in  this  scheme.  The  spatial  finite-difference  operator  is  the 
DRP/opt  MacConnack  scheme  (Hixon,  1997;  Zhang  and  Chen,  2006),  whose  formal 
accuracy  is  4th  order  with  optimization  in  numerical  dispersion  and  dissipation  errors;  the 
time  marching  scheme  can  be  2nd-order  Euler  or  4th-order  Runge-Kutta  based  on  the 
trade-off  between  accuracy  and  computational  efficiency.  This  implementation  enables  us 
to  accurately  simulate  seismic  waves  in  a  spherical  Earth  with  the  fluid  layer  and  3D 
variations  in  wave  speeds  and  attenuation.  A  number  of  numerical  experiments  have  been 
conducted  to  validate  the  spherical  finite-difference  code  by  comparing  the  simulation 
results  with  the  synthetics  calculated  by  normal-mode  summation  (Figure  2). 

At  local  to  near  regional  scales,  scattering  of  short-period  waves  by  topography 
becomes  important.  To  account  for  topography  effects,  we  have  implemented  a  3D 
version  of  the  traction  imaging  method  (Zhang  and  Chen,  2006)  to  simulate  wave 
propagation  in  the  presence  of  irregular  free  surface  (Zhang  et  al.,  2008).  The  spatial 
finite-difference  operator  is  the  DRP/opt  MacConnack  scheme  (Hixon,  1997);  the  time 
marching  scheme  is  4th-order  Runge-Kutta.  In  this  algorithm,  boundary-conforming 
grids  follow  surface  topography,  avoiding  the  artifacts  associated  with  the  staircase 
approximation  to  arbitrarily  irregular  surface.  For  example,  the  simulation  of  strong 
ground  motion  after  the  May  12,  2008  Wenchuan  earthquake  at  the  east  margin  of  the 
Tibetan  plateau  incorporates  a  3D  earth  model  (CRUST  2.0,  Bassin  et  al,  2000),  finite- 
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fault  rupture,  and  realistic  surface  topography  (Zhang  et  ah,  2008).  The  model  has  a 
dimension  of  700  km  x  300  km  x  80  km  in  the  fault-parallel  (N50°E),  fault-normal 
(N40°W),  and  vertical  directions,  respectively.  It  is  discretized  with  a  horizontal-grid 
spacing  of  500  m  and  variable  vertical-grid  spacing  from  200  m  near  the  surface  to  1.5 
km  near  80  km  depth.  The  simulation  provided  reliable  synthetic  waveforms  up  to  -0.33 
Hz  in  the  sediment  basin  and  up  to  -0.8  Hz  elsewhere.  One  of  the  conclusions  of  the 
simulation  is  that  rough  topography  above  the  rupture  fault  causes  wave  scattering, 
resulting  in  significantly  larger  peak  ground  motion  on  the  apex  of  topographic  relief  than 
in  the  valley.  Topography  and  scattering  also  reduce  the  wave  energy  in  the  forward 
direction  of  fault  rupture  but  increase  the  peak  ground  motion  in  other  parts  of  the  basin. 
These  results  suggest  the  need  for  a  localized  hazard  assessment  in  places  of  rough 
topography  that  takes  topographic  effects  into  account.  From  a  technical  point  of  view, 
this  simulation  was  carried  out  on  a  104  cpu-core  cluster,  demonstrating  that  it  is  feasible 
to  simulate  broadband  wave  propagation  at  regional  distances  on  a  moderate 
supercomputer. 

Finally,  we  note  that  the  perfectly  matched  layer  (PMF)  has  been  widely  used  in 
seismic  wave  modeling  to  absorb  waves  at  nonphysical  boundaries  surrounding  the 
computational  domain.  However,  the  traditional  PMF  performs  poorly  in  absorbing  low- 
frequency  waves,  evanescent  waves  or  waves  at  grazing  incidence.  We  developed  an 
unsplit-field,  complex  frequency-shifted  perfectly  matched  layer  implementation  using 
auxiliary  differential  equations  (CFS  ADE-PMF)  in  seismic  wave  modeling  (Zhang  and 
Shen,  2009).  One  of  the  advantages  of  the  CFS  ADE-PMF  implementation  is  that  it  is 
straightforward  to  use  with  a  higher-order  time  matching  scheme  to  achieve  a  higher- 
order  accuracy. 

3.3.  Source  Characterization  and  Construction  of  Hierarchical  FD-SGT 
Database 

Most  earthquake  source  studies  to  date  utilize  ID  earth  models  in  the  calculation  of 
synthetic  waveforms  or  Green’s  functions.  To  accommodate  3D  structural 
heterogeneities,  segments  of  body  and  surface  waves  are  time-shifted  before  waveform 
fitting  (e.g.,  Zhao  &  Helmberger,  1994;  Zhu  &  Helmberger,  1996).  In  places  with  highly 
heterogeneous  structure,  this  “cut-and-paste”  of  waveform  segments  may  not  fully 
account  for  the  complexity  of  wave  phenomena  caused  by  3D  structures  (e.g.,  wave 
focus ing/defocusing).  Furthermore  one  of  the  major  sources  of  errors  in  moment  tensor 
inversion  is  phase  skipping  between  the  observed  and  synthetic  waveforms.  When 
traveltime  delays  due  to  3D  velocity  heterogeneities  are  comparable  to  or  larger  than  the 
wave  period,  phase  skipping  becomes  a  challenging  problem. 

Advance  in  computation  has  made  the  use  of  3D  Green’s  function  in  earthquake 
source  inversion  an  increasingly  attractive  option  in  routine  operation.  Fiu  et  al.  (2004) 
determined  the  source  mechanisms  of  small  to  moderate  earthquakes  in  southern 
California  using  synthetic  waveforms  calculated  with  the  spectral-element  method. 
Because  the  derivatives  of  the  source  parameters  are  determined  numerically  by 
differentiating  synthetics  with  respect  to  the  source  parameters,  up  to  10  forward 
simulations  are  needed  for  each  earthquake  (6  moment  tensor  components,  latitude, 
longitude,  depth,  and  the  reference  location).  When  the  reference  location  is  far  from  the 
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true  location,  the  derivatives  of  the  location  parameters  may  not  adequately  reflect  the 
non-linear  variation  of  waveforms  as  a  function  of  location.  Taking  advantage  of  the 
source-receiver  reciprocity,  Zhao  et  al.  (2006)  introduced  the  use  of  strain  Green’s  tensor 
in  source  inversion  based  on  3D  reference  models.  We  extend  the  SGT-database 
approach  to  regional  and  teleseismic  scales,  so  we  are  able  to  use  any  combination  of 
local  and  teleseismic  waveforms  in  source  inversion  (Shen  et  al.,  2009).  We  fit  filtered 
and  time-shifted  wavefonns  directly.  For  a  given  source  location,  the  waveforms  are 
linear  functions  of  the  moment  tensor  elements,  which  can  be  solved  as  a  linear  inverse 
problem.  To  account  for  perturbations  in  earthquake  source  location,  we  carry  out  a  grid 
search  in  the  vicinity  of  the  reference  location. 

Figure  3  shows  the  selected  earthquakes  in  the  southeast  Tibetan  plateau  and  the 
stations  used  in  moment  tensor  inversion  (Shen  et  al.,  2009).  Wave  simulations  are 
carried  out  at  three  scales  (local/near  regional,  far  regional,  and  teleseismic),  resulting  in 
three  FD-SGT  databases,  one  for  each  wave  simulation  level.  An  updated  crust  and  upper 
mantle  model  (CUB  2.0,  Shapiro  &  Ritzwoller,  2002)  is  used  as  the  3D  reference  model. 
The  FD-SGT  moment  tensors  agree  well  the  global  centroid  moment  tensor  (GCMT) 
solutions  (Fig.  4),  with  one  notable  exception  at  the  eastern  Himalayan  syntaxis  (event 
03.230  in  Figure  4).  The  GCMT  solution  shows  a  strike-slip  event  at  33  km  depth,  while 
our  FDSGT  solution  is  a  shallow  (2±16  km)  thrust  event.  Figure  5  shows  an  example  of 
the  fit  between  the  observed  waveforms  and  the  synthetics  calculated  from  the  FD-SGT 
databases.  The  regional  stations,  the  3D  model  and  shorter-period  body  waves  than  in  the 
GCMT  solutions  tend  to  provide  tighter  constraints  on  the  source  location,  particularly 
the  depth.  The  lowest  magnitude  of  the  earthquakes  in  Figure  4  is  Mw  4.35,  though  the 
method  is  applicable  to  smaller  earthquakes.  The  SGT  approach  is  most  effective  when 
the  model  predicts  significantly  better  waveform  fit  than  ID  models.  As  the  crust  and 
mantle  velocity  model  improves,  one  may  extend  the  useful  surface  waves  to  shorter 
periods  and  be  able  to  detennine  the  moment  tensors  of  lower  magnitude  earthquakes. 

3.4.  Sensitivity  Kernels  Based  on  3D  Reference  Model  and  Full  Wave 
Tomography 

In  recent  years,  seismologists  have  started  to  replace  the  great-circle  path  of  a  surface 
wave  with  sensitivity  kernels  having  a  finite  width  and  symmetrically  distributed  across 
the  great-circle  path  (e.g.,  Yoshizawa  and  Kennett,  2002;  Ritzwoller  et  al.,  2002a)  or  a 
more  physically  realistic  kernel  calculated  using  ID  (e.g.,  Li  and  Romanowicz,  1996; 
Zhou  et  al.,  2004,  2006)  or  3D  (e.g.,  Liu  and  Tromp,  2006;  Zhang  et  al.,  2007)  reference 
models.  Similarly,  several  studies  (e.g.,  Hung  et  al.,  2004;  Chen  et  al.,  2007b;  Sigloch  et 
al.,  2008)  have  replaced  body-wave  ray  paths  with  “banana-doughnut”  sensitivities 
calculated  in  ID  (e.g.,  Dahlen  et  al.,  2000)  or  3D  (e.g.,  Tromp  et  al.,  2005;  Zhao  et  al., 
2005;  Liu  and  Tromp,  2006;  Zhang  et  al.,  2007)  models. 

To  date  the  inversions  of  P-  and  S-wave  speeds  have  been  carried  out  separately 
under  the  assumption  that  P  traveltimes  are  affected  only  by  the  P-wave  speed  of  the 
elastic  media  and  S  traveltimes  by  the  S-wave  speed.  The  practice  may  stem  from  ray 
theory,  in  which  P  and  S  waves  are  completely  decoupled  except  for  conversions  at 
discontinuities.  Based  upon  numerical  and  analytical  solutions,  we  have  shown  that  for 
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finite-frequency  waves,  S-speed  perturbation  may  have  significant  effects  on  P 
waveforms  (Zhang  and  Shen,  2008).  At  places  near  the  source  or  receiver,  the  cross¬ 
dependence  of  P  waveforms  on  S-speed  perturbation  is  substantial  (Figure  6).  For 
example,  the  integrated  absolute  P-traveltime  sensitivity  to  S-speed  perturbation  is  31% 
of  that  to  P-speed  perturbation  when  the  source-receiver  distance  is  seven  times  the 
wavelength  (LA=7,  where  X  is  the  wavelength  and  L  the  distance  to  the  receiver).  At  a 
source-receiver  distance  of  100  wavelengths,  the  ratio  is  about  5%.  This  finite-frequency 
effect  has  important  implications  for  seismic  tomography,  particularly  at  regional  and 
local  scales,  where  target  structures  are  often  located  within  several  tens  of  wavelengths 
from  receivers  or  from  both  sources  and  receivers.  Even  for  regions  far  from  sources  and 
receivers,  neglecting  this  cross-dependence  means  that  the  estimates  of  the  P  speed  near 
the  source  and  receiver  are  likely  biased  and  such  errors  may  propagate  into  the  rest  of 
the  model  in  tomographic  inversion. 

Figure  7  shows  the  phase  sensitivities  of  a  regional  Rayleigh  wave  to  S-  and  P-wave 
speeds.  Similar  to  the  results  for  southern  California  in  Liu  and  Tromp  (2006),  there  are 
large  3D  variations  in  the  finite-frequency  kernels  due  to  wave-speed  variations. 
Although  commonly  neglected  in  surface  wave  tomography,  P-wave  speed  affects 
Rayleigh  waves  (also  see  Zhou  et  al.  (2004)  and  Liu  and  Tromp  (2006)).  The  maximum 
sensitivity  to  P-wave  speed  is  on  the  surface  and  at  a  level  comparable  to  the  Vs  kernel. 
It  decreases  to  -10%  of  the  maximum  at  depth  of  -Z/4.  Since  the  peak  of  the  Vs 
sensitivity  is  at  the  depth  of  -Z/3  and  heterogeneities  in  the  mantle  are  generally  much 
weaker  than  in  the  crust,  it  may  be  necessary  to  include  the  sensitivities  to  P-speed  and 
density  for  long-period  (-50+  s)  Rayleigh  waves  in  a  high-resolution  inversion.  When 
integrated  with  P/Pn,  Rayleigh  waves  may  provide  additional  constraints  on  the  P-wave 
structure  in  the  crust.  At  the  very  least,  incorporating  the  effects  of  P-speed  on  Rayleigh 
waves  may  minimize  a  likely  cause  contributing  to  the  apparent  Love-Rayleigh 
discrepancy  and  radial  anisotropy. 

Figure  8  shows  the  P  wave  velocity  model  of  the  southeast  Tibetan  plateau  in  the 
shallow  (8  km)  crust.  The  model  is  derived  from  a  joint  inversion  of  -10,000 
measurements  of  phase  anomalies  of  P/Pn,  S/Sn  and  Rayleigh/Love  waves  in  three 
frequency  bands  ranging  from  0.008-0.1  Hz.  After  two  iterations,  the  combined  velocity 
anomalies  relative  to  the  initial  reference  model  (CUB  2.0,  Shapiro  and  Ritzwoller,  2002) 
show  that  the  observations  require  a  band  of  velocity  increase  just  north  of  the  Indus  Yalu 
suture.  Another  band  of  high  velocity  anomalies  is  located  at  longitude  ~102°E.  These 
high-velocity  anomalies  correlate  well  with  areas  of  granitic  and/or  more  metamorphosed 
Proterozoic  rocks  among  the  regional  Mesozoic  background  (Geological  map  of  the 
Tibetan  plateau,  China  Academy  of  Geological  Sciences,  Chengdu).  Compared  to 
previous  models  and  tomography  of  similar  scales,  the  unprecedented  P-wave  resolution 
at  such  a  shallow  depth  is  attributed  in  part  to  the  sensitivity  of  Rayleigh  waves  to  P  wave 
speed  perturbations. 

3.  DISCUSSION  AND  CONCLUSIONS 

We  have  systematically  developed  a  finite-difference  method  to  accurately  and 
efficiently  simulate  seismic  wave  propagation  at  local  to  regional  scales.  The  non- 
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staggered-grid  finite-difference  code  in  the  spherical  coordinate  system  is  well  suited  for 
regional  and  global  wave  propagation,  while  the  finite-difference  method  with  boundary- 
conforming  grids  that  follow  topography  provides  higher  accuracy  at  local  and  near 
regional  scales.  The  CFS  ADE-PML  implementation  of  the  Perfectly  Matched  Layer 
makes  it  possible  to  use  a  relatively  small  number  of  boundary  layers  to  efficiently  absorb 
waves  incident  at  a  near-grazing  angle,  low  frequency  waves  and  evanescent  waves. 

The  local-  and  global-scale  finite-difference  methods  make  it  possible  to  construct 
hierarchical  FD-SGT  databases  with  great  flexibility.  These  SGT  databases  can  be  used 
to  invert  for  the  source  moment  tensors  and  locations,  a  highly  nonlinear  problem,  in  a 
global-minimum-misfit  search.  Because  the  SGT  databases  are  pre-calculated,  the 
inversion  of  the  source  parameters  does  not  involve  forward  wave  simulations  “on  the 
fly”  and  thus  can  be  carried  out  in  real  time  or  near  real  time. 

The  full-wave  tomography  method  provides  an  accurate  way  to  integrate  body-  and 
surface-wave  observations  in  a  self-consistent  way  via  3D  sensitivity  kernels.  Our  ability 
to  update  the  3D  models  means  we  can  linearize  the  tomographic  problem  by 
progressively  detennining  heterogeneities  at  various  scales  from  broad  and  smooth 
features  to  small-scale,  large  and  sharp  variations.  For  earthquake  tomography,  the  rate 
of  convergence  depends  on  the  initial  reference  model,  since  it  affects  the  source  moment 
tensor  solutions,  and  thus  measurements  of  the  differences  between  the  observed  and 
synthetic  waveforms  (e.g.,  traveltime).  One  possible  strategy  to  deal  with  this  nonlinear 
problem  is  to  use  Empirical  Green’s  Functions  (EGFs)  extracted  from  ambient  noise  to 
obtain  a  robust  model  for  the  crust  and  shallow  upper  mantle  before  the  inversion  for  the 
earthquake  source  moment  tensors  and  subsequent  earthquake  tomography.  Because 
EGFs  are  dominated  by  fundamental-mode  Rayleigh  and  Love  waves,  they  can  be  used 
in  the  same  way  in  full-wave  tomography  as  earthquake  Rayleigh  and  Love  waves. 

4.  RECOMMENDATIONS 

An  increasing  number  of  seismologists  have  appreciated  the  advantages  of  the  full- 
wave  approach  based  on  3D  reference  models  although  only  a  few  groups  have  actually 
used  it  to  date.  This  is  partly  due  to  the  fact  that  the  method  requires  the  use  of  high- 
performance  supercomputers.  One  promising  area  of  developments  in  computational 
technology  is  the  use  of  clusters  constructed  with  Graphic  Processing  Units  (GPU), 
which  can  improve  finite-difference  simulations  by  a  factor  of  10  to  30  compared  to 
CPU-based  clusters  (Abdelkhalek  et  al.,  2009).  So  a  well-adapted  computation  code  on  a 
moderate  ($100K)  GPU  cluster  can  be  as  powerful  as  a  million-dollar  CPU  cluster.  The 
porting  of  finite-difference  codes  to  GPU  clusters  should  be  relatively  straightforward 
and  inexpensive.  This  will  significantly  lower  the  costs  of  computational  resources 
needed  to  do  full-wave  source  moment  tensor  and  tomography  inversion. 
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Figure  2  Schematic  illustration  of  the  unified  tomography  and  moment  tensor  inversions  based 
on  3D  FD-SGT  databases. 
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Figure  3  Comparisons  between  the  3-component  synthetics  calculated  by  the  spherical  finite- 
difference  algorithm  (dashed  line)  and  normal-mode  summation  (solid  line).  The  double  couple 
source  with  the  strike,  dip  and  rake  of  45°,  45°  and  90°,  respectively,  is  located  at  (0°N,  3°E).  The 
receiver  is  at  (0°N,  18°E)  on  the  surface,  15°  east  of  the  source.  The  fit  of  the  two  solutions  is 
excellent,  verifying  the  accuracy  of  the  finite-difference  calculation. 
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Figure  4  Selected  earthquakes  (dots)  in  the  southeast  Tibetan  plateau  during  2001-2004  and  the 
stations  (triangles)  used  in  moment  tensor  inversions.  White  rectangles  outline  the  three 
geographic  scales  of  finite-difference  wave  simulation  (1)  local/near  regional,  (2)  far-regional, 
and  (3)  teleseismic.  Shen  et  al.  (2009). 
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Figure  5  Comparison  of  the  FD-SGT  (white-and-blue  symbols)  and  GCMT  (white-and-red 
symbols)  solutions  for  13  earthquakes  during  2001-2004  in  the  southeast  Tibetan  plateau.  The 
numbers  above  the  FDSGT  solution  denote  the  year  and  Julian  day  of  the  earthquake.  Also 
plotted  are  the  source  mechanisms  of  10  additional  earthquakes  without  GCMT  solutions  (black- 
and-white  symbols). 
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Figure  6  (a-d)  Observed  waveforms  (thick  gray  lines)  for  event  01.300  of  Figure  4  are 
compared  with  the  synthetics  for  the  best  source  moment  tensor  solution  (red  lines).  The 
frequency  bands  are,  from  left  to  right,  0.00833-0.0167,  0.0167-0.0333,  0.0333-0.0667, 
and  0.0667-0.167  Hz.  Dark  thick  lines  are  time-shifted  observed  waveforms  used  in  the 
final  inversion.  For  the  first  arrivals  (in  c  and  d),  the  time  is  reduced  to  align  the  arrivals 
within  the  time  window.  The  amplitude  is  normalized  by  individual  trace,  (e) 
Normalized  waveform  misfit  as  a  function  of  depth.  Crosses  mark  the  grid-search 
depths.  The  shaded  region  indicates  the  1-a  confidence  limit,  (f)  Contours  of  various 
confidence  limits  at  the  depth  of  the  global  best  solution  (8  km  for  this  event).  The  cross 
marks  the  FDSGT  location,  the  square  the  EHB  location  (E.R.  Engdahl,  per.  comm.),  and 
the  diamond  the  global  centroid  moment  tensor  location.  From  the  innermost  contour 
surrounding  the  FDSGT  solution,  the  contours  are  the  68%  (dark  line),  80%,  90%,  95% 
(dark  line),  and  99%  confidence  limits,  (g)  Triangles  and  cross  mark  the  stations  used  in 
the  source  inversion  and  the  event,  respectively.  Gray  contours  follow  0,  2000  and  4000 
m  elevation.  Shen  et  al.  (2009). 
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Figure  7  The  relative  magnitudes  of  the  integrated  absolute  P  traveltime  (asterisk)  and  amplitude 
(triangle)  sensitivities  to  S-  and  P-speed  perturbations  vary  with  the  normalized  distance  between 
the  source  and  receiver.  From  Zhang  and  Shen  (2008). 
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Figure  8  Phase  sensitivity  kernels  of  a  regional  Rayleigh  wave  to  S  wave-speed  (Top)  and  P 
wave-speed  (Bottom)  perturbations  relative  to  a  3D  reference  model.  The  period  of  the  wave  is 
10-20  s.  The  two  vertical  cross-sections  cut  the  model  along  the  great  circle  path  between  the 
earthquake  source  (in  the  southeast  comer  of  the  model)  and  receiver  (in  the  northwest  comer  of 
the  model)  and  perpendicular  to  the  great  circle.  Also  shown  is  a  quadrant  of  a  horizontal  section 
on  the  surface.  The  scale  is  ±2xl0'13  s/m3  for  both  S  and  P  wave  speed. 
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Figure  9  P  wave  velocity  variations  at  8  km  depth  in  the  southeast  Tibetan  plateau.  The  initial 
3D  reference  model  (top  right)  is  converted  from  CUB  2.0  by  a  Vp/Vs  ratio  of  1.74  in  the  crust. 
The  unit  of  the  color  bar  is  km/s.  The  left  panels  show  the  velocity  perturbations  after  the  1st  and 
2nd  iteration  and  the  total  velocity  change  in  percent.  The  middle-right  and  bottom-right  panels 
shows  the  updated  models  after  the  1st  and  2nd  iterations.  The  high-velocity  anomalies  north  of 
the  Indus-Yalu  suture  (IYS)  and  Bangong-Nujiang  suture  (BNS)  and  along  the  Xianshuihei  Fault 
system  (XF)  are  correlated  with  areas  of  granitic  and/or  older  and  more  metamorphosed 
Proterozoic  rocks. 
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Running  head:  ADE  CFS-PML  for  Seismic  Wave  Modeling 

ABSTRACT 

The  complex  frequency-shifted  PML  (CFS-PML)  is  more  efficient  than  the  standard  PML  to  absorb 
low-frequency  waves,  evanescent  waves  and  near-grazing  incident  waves.  In  seismic  wave  model¬ 
ing,  the  CFS-PML  has  been  implemented  by  the  first-order-accuracy  convolutional  PML  (CPML) 
technique  or  second-order-accuracy  recursive  convolution  PML  (RIPML)  technique.  Both  involve 
an  algorithm  different  from  the  numerical  scheme  used  in  the  interior  domain  to  update  the  auxil¬ 
iary  memory  variables  and  thus  can  not  be  directly  used  in  higher-order  time  marching  schemes. 
In  this  paper,  we  introduce  an  alternative  unsplit-field  CFS-PML  implementation  using  auxiliary 
differential  equations  (ADE)  to  update  the  auxiliary  memory  variables.  This  ADE  CFS-PML  im- 
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plementation  results  in  complete  first-order  differential  equations  for  both  the  wavefield  components 
and  the  auxiliary  memory  variables.  Thus  the  same  numerical  scheme  used  in  the  interior  domain 
can  be  used  to  solve  the  ADE  CFS-PML  equations.  We  implement  the  ADE  CFS-PML  in  both 
the  FDTD  method  and  a  nonstaggered-grid  MacCormack  DRP/opt  finite-difference  method  with 
the  4th-order  Runge-Kutta  scheme  to  demonstrate  its  straightforward  implementation  in  different 
numerical  time  marching  schemes.  We  perform  a  series  numerical  tests  to  investigate  the  optimal 
values  of  the  CFS-PML  parameters  in  seismic  wave  modeling  and  show  that  the  scaling  factor  /? 
transforms  the  PML  to  a  transversely  isotropic  material.  For  normal  incident  waves,  (5  only  delays 
wave  propagation  but  does  not  increase  damping.  For  oblique  incident  waves,  it  bends  the  wave- 
front  towards  the  normal  direction  of  the  PML  layer,  thus  increases  damping.  Theoretically,  a  larger 
/?  can  bend  the  wavefront  more  efficiently,  but  a  very  large  (3  may  cause  strong  spurious  reflec¬ 
tions  if  it  causes  the  wavelength  to  be  much  shorter  than  what  the  numerical  scheme  can  resolve. 
Our  numerical  tests  indicate  that  the  optimal  f3  reduces  the  points  per  dominant  wavelength  at  the 
outermost  boundary  to  three,  around  half  the  value  required  by  the  numerical  scheme. 

At  the  free  surface,  the  PML  equations  should  be  derived  taking  the  free  surface  boundary 
condition  into  account  in  finite-difference  methods.  Otherwise,  the  free  surface  in  the  PML  Layer 
will  cause  instability  in  the  nonstaggered-grid  finite-difference  scheme  or  ineffective  absorption  of 
surface  waves  when  the  free  surface  coincides  with  the  normal  stress  components  in  the  staggered- 
grid  finite-difference  scheme. 

Our  numerical  test  of  a  two-layered  medium  with  surface  topography  in  a  narrow  vertical  slice 
indicates  that  we  can  use  a  narrow  slice  mesh  with  the  ADE  CFS-PML  to  efficiently  simulate  full 
wave  propagation  in  complex  structure  models. 
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INTRODUCTION 

Numerical  study  of  seismic  wave  propagation  usually  needs  to  simulate  waves  propagating  in  an 
unbounded  space.  However,  numerical  simulation  is  often  implemented  in  a  finite  computational 
domain  in  practice.  Thus  some  special  treatments  are  needed  at  the  computational  boundaries  to 
absorb  waves  propagating  outward.  Techniques  based  on  two  concepts  are  commonly  used,  one 
is  called  Absorbing  Boundary  Conditions  (ABCs),  which  impose  a  proper  boundary  condition  at 
the  outermost  boundaries  to  satisfy  the  condition  that  waves  only  propagate  outward  (Clayton  and 
Engquist,  1977;  Higdon,  1986,  1990;  Randall,  1988;  Liao  et  al.,  1984;  Bayliss  et  al.,  1986),  the 
other  can  be  grouped  into  Absorbing  Boundary  Layers  (ABLs,  or  sponge  layers),  which  use  finite 
layers  to  gradually  damp  wave  amplitude,  so  using  the  Dirichlet  boundary  condition  at  the  outermost 
boundary  does  not  generates  a  strong  reflection  (Cerjan  et  al.,  1985;  Sochacki  et  al.,  1987).  Strengths 
and  shortcomings  of  each  group  can  be  found  in  Lesta  and  Vilotte  (2005)  and  Konratitsch  and  Martin 
(2007).  Perfectly  matched  layer  (PML)  (Berenger,  1994),  one  of  the  ABLs  techniques,  is  proven  to 
effectively  absorb  waves  at  a  wide  range  of  incident  angles  with  only  several  to  tens  layers  and  has 
become  widely  used  in  elastic  wave  modeling  (Chew  and  Liu,  1996;  Hastings  et  al.,  1996;  Collino 
and  Tsogka,  2001;  Marcinkovich  and  Olsen,  2003;  Wang  and  Tang,  2003). 

The  original  PML  in  Berenger  (1994)  (we  refer  to  this  version  as  the  standard  PML  in  this  pa¬ 
per)  adopted  a  non-physical  splitting  of  the  wave  field  components  and  equations,  which  leads  to 
two  different  sets  of  equations  for  the  interior  domain  and  the  PML  layer,  and  thus  are  not  trivial 
to  incorporate  into  an  existing  numerical  modeling  code.  Lurthermore,  the  split-field  PML  is  math¬ 
ematically  weakly  well-posed  (Abarbanel  and  Gottlieb,  1997).  Alternative  intepretations  of  the 
PML  as  an  artificial  anisotropic  medium  (Sacks  et  al.,  1995;  Gedney,  1996),  or  complex  coordinate 
stretching  (Chew  and  Weedon,  1994;  Teixeira  and  Chew,  2000)  led  to  unsplit-field  PML  implemen- 
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tations  involving  convolution  terms  (Wang  and  Tang,  2003;  Komatitsch  and  Martin,  2007),  integral 
terms  (Zeng  and  Liu,  2004;  Drossaert  and  Giannopoulos,  2007b),  or  auxiliary  differential  equations 
(Ramadan,  2003).  Hagstrom  (2003)  introduced  a  modified  modal  solution  technique  to  derive  PML 
equations  and  was  used  by  Appelo  and  Kreiss  (2006)  in  2D  elastic  wave  modeling.  This  technique 
can  provide  a  stable  PML  implementation  by  choosing  the  parameters  satisfying  the  geometric  rela¬ 
tions  in  anisotropic  media  (Becache  et  al.,  2003).  Recently,  Meza-Fajardo  and  Papageorgiou  (2008) 
showed  that  stable  result  can  be  obtained  in  anisotropic  media  using  a  3D  damping  profile  in  the 
PML  based  on  the  concept  of  complex  coordinate  stretching. 

Though  the  standard  PML  has  been  widely  used  in  seismic  wave  modeling,  it  can  generate  large 
spurious  reflections  for  near-grazing  incident  waves,  low  frequency  waves  or  evanescent  waves 
(Festa  and  Vilotte,  2005;  Komatitsch  and  Martin,  2007;  Drossaert  and  Giannopoulos,  2007a).  The 
complex  frequency  shifted  PML  (CFS-PML)  (Kuzuoglu  and  Mittra,  1996)  is  more  efficient  in  such 
circumstances  using  a  frequency  dependent  damping  (Festa  and  Vilotte,  2005;  Festa  et  al.,  2005;  Ko¬ 
matitsch  and  Martin,  2007;  Drossaert  and  Giannopoulos,  2007b, a).  The  CFS-PML  was  originally 
implemented  in  a  split-field  form  and  three  auxiliary  memory  variables  are  needed  for  each  deriva¬ 
tive  (Gedney,  1998).  Roden  and  Gedney  (2000)  derived  an  unsplit-field  CFS  Convolutional-PML 
(C-PML)  equation  involving  a  convolution  term  that  can  be  efficiently  calculated  using  a  recursive 
convolution  algorithm  (Luebbers  and  Flunsberger,  1992).  This  technique  has  been  used  in  seis¬ 
mic  wave  modeling  in  elastic  media  (Komatitsch  and  Martin,  2007;  Drossaert  and  Giannopoulos, 
2007a)  and  poroelastic  media  (Martin  et  al.,  2008a).  Drossaert  and  Giannopoulos  (2007b)  derived 
an  alternative  CFS-PML  implementation  that  involves  integral  terms,  which  can  be  calculated  by  a 
recursive  integration  algorithm.  The  recursive  integration  in  C-PML  is  of  first-order  accuracy  (Gi¬ 
annopoulos,  2008),  and  the  trapezoidal  rule  in  RIPML  is  of  second-order  accuracy  (Drossaert  and 
Giannopoulos,  2007b).  Both  approaches  need  an  algorithm  different  from  the  numerical  scheme 
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used  in  the  interior  domain  to  update  the  auxiliary  memory  variables  and  thus  can  not  be  directly 
used  with  higher-order  time  marching  schemes,  such  as  a  4th-order  Runge-Kutta  scheme  used  in 
the  nonstaggered  finite-difference  velocity-stress  method  (Zhang  and  Chen,  2006)  on  a  boundary- 
conforming  grid  to  simulate  seismic  wave  propagation  in  the  presence  of  surface  topography,  in 
which  the  higer-order  multi-stage  Runge-Kutta  scheme  is  important  to  allow  a  larger  time  step 
(Hixon,  1997).  In  the  field  of  electromagnetic  simulation,  Ramadan  (2003)  proposed  an  unsplit- 
field  implementation  of  the  standard  PML  using  auxiliary  differential  equations  (ADE),  and  Wang 
and  Liang  (2006)  extended  this  approach  to  the  CFS-PML  with  a  2D  alternating  direction  implicite 
(ADI)  finite-difference  time-domain  (FDTD)  method.  This  ADE  CFS-PML  implementation  results 
in  complete  first-order  differential  equations  for  both  the  wavefield  components  and  the  auxiliary 
memory  variables,  thus  the  same  numerical  scheme  used  in  the  interior  domain  can  be  used  to  solve 
the  ADE  CFS-PML  equations  no  matter  whether  a  2nd-order  leap-frog  scheme  or  a  higher-order 
Runge-Kutta  scheme  is  used. 

In  this  paper,  we  implement  the  ADE  CFS-PML  in  the  seismic  wave  simulation  both  in  the 
FDTD  method  and  the  nonstaggered-grid  finite-difference  method  (Zhang  and  Chen,  2006)  using 
a  4th-order  Runge-Kutta  time  marching  scheme.  We  derive  the  ADE  CFS-PML  equations  for  the 
velocity-stress  equations.  Then  we  analyze  the  free  surface  boundary  condition  inside  the  PML 
layer  and  modify  the  PML  equations  based  on  the  free  surface  condition  to  avoid  the  instability 
caused  by  the  incompatibility  between  the  PML  equations  and  the  free  surface  boundary  condition. 
We  analyze  the  role  of  the  scaling  factor  (3  in  the  CFS  stretching  function,  which  is  still  not  clear  in 
the  literature  of  seismic  wave  modeling.  We  perform  a  series  of  numerical  tests  in  a  thin  slab  model 
to  validate  the  ADE  CFS-PML  implementation  and  investigate  the  optimal  values  of  the  parameters 
in  the  CFS  stretching  function.  We  add  the  free  surface  to  the  thin  slab  model  to  demonstrate  the 
efficiency  of  the  modified  PML  equations  taking  into  account  the  free  surface  condition  in  absorbing 
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surface  waves.  Finally,  we  simulate  seismic  wave  propagation  in  a  two-layered  medium  with  surface 
topography  in  a  narrow  vertical  slice  to  demonstrate  the  implementation  of  the  ADE  CFS-PML  in 
the  4th-order  Runge-Kutta  scheme  and  its  potential  effectiveness  in  seismic  exploration. 

VELOCITY-STRESS  EQUATION  AND  FINITE  DIFFERENCE  NUMERICAL 

SCHEME 

We  consider  the  following  velocity-stress  formulation  of  the  elastic  wave  equations  for  an  isotropic 
medium  in  the  Cartesian  coordinate  (x,  y.  z), 


v,i  =  -V  ■  a, 

(1) 

P 

aj  =  c:  [Vv  +  (Vv)T]  , 

(2) 

where  v  =  (vx,vy,vz)T  is  the  velocity  vector,  a  is  the  stress  tensor  and  the  components  are 
atJ,  i.  j  £  (x,  y.  z"),  p  is  the  mass  density,  and  c  is  the  stiffness  tensor.  In  this  paper,  a  comma 
followed  by  t,  x,  y  and  2  in  subscript  means  a  derivative  with  respect  to  t,  x,  y  or  2.  Equation  1-2 
can  also  be  written  in  a  single  vector  form  appropriate  for  the  4th-order  Runge-Kutta  scheme, 

W  >t  =  AW,,  +  BW,j  +  CW,2)  (3) 

where  W  =  (vx,  vy .  vz ,  axx ,  ayy,  azz,  oyz ,  axz .  axy )  is  the  velocity-stress  vector.  A,  B  and  C  are 
the  coefficient  matrices  that  can  be  easily  derived.  Sometimes  it  is  more  convinenient  to  discuss  the 
PML  implementation  in  the  component  form,  thus  we  express  the  vx  component  of  equation  1 , 

PVx,t  —  &xx,x  ~P  &xy,y  "F  &xz,z-  (4) 

We  use  two  numerical  schemes  to  solve  the  above  velocity-stress  equations,  one  is  the  finite- 
difference  time-domain  (FDTD)  method  with  a  2nd-order  accuracy  leap-frog  time  marching  scheme 
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on  a  staggered  grid  (Levander,  1988;  Graves,  1996),  the  other  is  the  nonstaggered-grid  finite- 
difference  method  with  a  4th-order  Runge-Kutta  time  marching  scheme  on  a  boundary-conforming 
grid  (Zhang  and  Chen,  2006).  Since  PML  implementations  do  not  rely  on  how  spatial  derivatives 
are  calculated,  we  only  list  the  time  marching  schemes  here.  The  2nd-order  leap-frog  scheme  can 
be  written  as 


■”+1/2  =  (T71-!/2  +  AfL(Vn), 

(5) 

vn+1  =  vn  +  AfL(crn+1/2), 

(6) 

where  superscript  n  —  1/2,  n,  n  +  1/2  and  n  +  1  means  time  level,  L(<7n+1/2)  is  the  right  hand  side 
of  the  equation  2  evaluated  at  time  level  n  +  1/2  using  an  arbitary  order  accuracy  spacial  staggered 
finite-difference  operator,  L(vn)  is  that  of  the  equation  1  at  time  level  n.  The  4th-order  Runge-Kutta 
scheme  can  be  expressed  as 

/i(1)  =  AtL(Wn)  (7) 

h{2)  =  AtL(Wn  +  a2h^) 
h{3)  =  AtL(Wn  +  a3h,W) 
h (4)  =  Ath{Wn  +  aAh^) 

Wn+1  =  Wn  +  ft/jW  +  /32h(2)  +  /33fi(3)  +  /34/i(4), 

where  a2  ~  0:4  and  3\  ~  /)4  are  the  coefficients  of  the  4th-order  Runge-Kutta  scheme,  L(ii/n) 
means  the  right  hand  side  of  equation  3  calculated  for  variable  Wn  using  the  spacial  operator.  Here 
the  DRP/opt  MacCormack  operator  (Hixon,  1997)  is  used. 
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COMPLEX  FREQUENCY-SHIFTED  PML  USING  AUXILIARY  DIFFERENTIAL 

EQUATIONS 


Each  face  of  a  3D  computational  grid  can  be  terminated  by  a  PML  layer  (Figure  1).  In  this  paper,  we 
will  discuss  the  PML  at  the  positive  x  face.  The  PML  at  the  opposite  face  and  along  other  directions 
can  be  derived  in  a  similar  way.  For  simplicity  of  notations,  we  also  suppose  the  PML  at  the  positive 
x  face  starts  from  x  =  0. 


Based  on  the  concept  of  complex  coordinate  stretching  (Chew  and  Weedon,  1994),  the  equations 
inside  the  PML  layer  have  exactly  the  same  form  as  that  (equations  1-2)  in  the  physical  domain 
except  that  the  coordinate  x  is  replaced  with  a  complex  stretched  coordinate  x, 

nx 

X=  sx(r])dr],  (8) 

■Jo 

where  sx  is  the  complex  stretching  function  that  determines  the  characteristics  of  the  PML  model. 
From  relationship, 


d  Id 
dx  sx  dx  ’ 


(9) 


one  can  transform  equations  in  the  complex  coordinate  (x,  y,  z)  into  the  Cartesian  coordinate  (x.  y.  z) 
where  the  numerical  scheme  will  be  used.  Using  equation  4  as  an  example,  we  have 


Ida  XX  ,  daxy  ,  da, 
1 0JpVx  = - ^ - 1 - 70^  + 


dx 


dy 


dz  ’ 


(10) 


where  a  hat  over  a  quantity  denotes  the  value  in  the  frequency  domain.  Once  sx  is  determined, 
we  can  transform  this  equation  into  time  domain  to  get  the  PML  equation  suitable  for  time  domain 
numerical  methods.  The  special  formulation  of  sx  is  the  key  to  distinguishing  the  standard  PML,  in 
which 


sx{x) 


dx(x) 

1  + 

it o 

8 


(11) 
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and  the  CFS-PML,  in  which 


/  X  a  ,  N  ,  4  (X) 

Sx{x)  =  (3x{x)  H -  .  . 

ax[x)  + l u 


(12) 


In  above  equations,  dx  >  0  is  the  attenuation  factor  that  causes  the  amplitude  of  the  wavefield  to  be 
exponentially  reduced  inside  the  PML  layer,  ax  >  0  is  the  frequency  shifted  factor  that  makes  the 
attenuation  frequency  dependent,  j3x  >  1  is  the  scaling  factor,  which  has  been  found  in  numerical 
tests  to  be  important  for  absorption  of  evanescent  waves  (Liu,  1999;  Drossaert  and  Giannopoulos, 
2007a)  and  near-grazing  incident  waves  (Drossaert  and  Giannopoulos,  2007b),  but  the  reason  has 
not  been  theoretically  explained  in  the  literature  of  seismic  wave  modeling,  we  will  discuss  the 
role  of  (3X  in  later  sections.  dx(x),  ax(x),  fix(x)  and  therefore  sx(x)  are  all  functions  of  x.  In  the 
following,  we  omit  (x)  for  brevity. 


The  basic  idea  of  the  ADE  implementation  of  the  CFS-PML  is  to  separate  the  derivative  with 
respect  to  the  complex  coordinate  into  two  parts  and  only  keep  ic o  term  in  the  denominator  in  one 
of  them. 


1 


1 


1  1 


dl 


fix  +  fix  fix  (ax  4-  ioj)f3x  +  dx 


Thus, 


1  O&xx  ^  1  dx 

sx  dx  (3X  °xx'x  px  (ax  +  i u)f3x  +  da 


Defining  an  auxiliary  memory  variable  T£x> 


dx 


TiX  - 

XX  —  /  ,  .  \  n  ,  7  VXx,X-) 

{cxx  T  iw)/3ac  T  dx 


we  can  write  equation  14  as, 


1  ddxx  _  1  „  1  rf, x 

SX  dx  -  p/XX’X  fixXX' 


(13) 


(14) 


(15) 


(16) 
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Taking  this  equation  into  equation  10,  we  obtain  the  unsplit-field  CFS-PML  equation  for  the  vx 
component  in  the  frequency  domain, 

i ujpvx  — ~7j-&xx,x  ~o~'^'xx  ®xy,y  ®xz,zi  (1*7) 

Px  Px 

where  the  memory  variable  Txx  is 

i uT%x  +  [ax  +  px  =  ^axx,x.  (18) 

Transforming  equations  17  and  18  to  the  time  domain,  we  get  the  final  ADE  CFS-PML  equation  of 

vx  as 

PVx,t  - ®XX,X  T  CTXy  y  +  <JXZ,Z  T 

and 

n.r.t  +  («X  +  T-;x  =  (20) 

It  is  worth  noting  that  we  write  the  time  domain  PML  equation  as  a  correction  to  the  original  wave 
equation  because  it  is  easy  to  achieve  parallel  balance  and  implement  in  an  existing  code  (Drossaert 
and  Giannopoulos,  2007a).  The  ADE  CFS-PML  equations  for  axis  y  and  z  and  other  components 
can  be  obtained  in  the  same  manner.  The  complete  ADE  CFS-PML  equations  for  the  velocity-stress 
equations  are  in  Appendix  A. 


<7 


rpx 

xx, x  a  J-xxi 

Px 


(19) 


Since  the  velocity-stress  equations  1-2  and  the  ADE  CFS-PML  equations  19-20  are  all  first- 


order  partial  differential  equations  with  respect  to  time  t,  they  can  be  updated  by  the  same  time¬ 


stepping  scheme.  In  the  staggered  2nd-order  leap-frog  scheme,  equation  20  can  be  dicretized  as 


r,x|n-|-l 


_  rpx\n 
1  XX 


At 


+ 


rpx\n+l  rpX \n  , 

J-xx  "i  xx  _  (*x  n+ 1/2 

~2  =  ~pxaxx’x 


(21) 


Then  we  can  update  Txx  and  vx  through 


rpx\n-\-l 
■*-  XX 


2  -  A t(ax  +  dx/f3x)  |n  2A tdx//3x 

_  n+l/2 

2  +  A t(ax -\- dx / (3X)  11  2  +  At(ax  +  dx/f3x) 


(22) 
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and 


n+1  _  n 


=<  + 


((Tn+1/2 
p  xx,x 


+  CT 


n+1/2 

xy,y 


a 


At 

P 


n+l/2\  ,  ^  f  J_  _  1 

XZ’Z  ’  P  \Px 

A  tdx/(3x 


a 


n+ 1/2 
XX, X 


2  +  A  t{ax  +  dx/ (3X) 


rpx\n 

XT. 


2  +  A  t(ax  +  dx/ (3X) 


<7 


n+1/2 


■ 


(23) 


The  implementation  in  the  4th-order  Runge-Kutta  scheme  is  also  straightforward.  One  can  just 
includes  the  memory  variables  in  the  solution  vector  W  in  equation  3,  then  uses  the  Runge-Kutta 
scheme  (equation  7)  to  update  the  PML  equations  including  the  memory  variable  equations. 


The  computer  memory  requirement  of  the  ADE  CFS-PML  is  the  same  as  the  C-PML  (Ko- 
matitsch  and  Martin,  2007)  or  RIPML  (Drossaert  and  Giannopoulos,  2007b),  one  memory  variable 
per  normal  derivative  to  the  PML  layer. 


ADE  CFS-PML  BOUNDARY  AND  FREE  SURFACE  BOUNDARY  CONDITION 

To  our  knowledge,  there  is  no  publcation  on  how  to  implement  the  free  surface  boundary  condition 
in  the  PML  region,  which  is  crucial  to  efficiently  absorb  surface  waves  in  the  FDTD  method  or 
avoid  the  stability  problem  in  the  nonstaggered  finite-difference  method. 

In  the  physical  region,  the  free  surface  boundary  condition  requires  the  traction  to  be  zero, 

n  •  a  =  0,  (24) 

where  n  is  the  outward  normal  vector  of  the  free  surface.  For  a  flat  surface,  the  zero-traction 
condition  is  reduced  to  a  zero-stress  condition, 

7  zz  =  0,  7yZ  =  0,  7Xz  =  0.  (“5) 

In  the  PML  layer,  the  zero- valued  stress  components  should  remain  zero  at  the  free  surface  since 
the  effect  of  the  PML  is  to  exponentially  damp  the  wavefield  to  zero. 
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Because  the  velocity  derivatives  are  related  to  the  stress  components  through  the  stress-strain 
relation  (equation  2),  the  zero-stress  condition  leads  to  velocity  derivative  constraints  at  the  free 
surface, 


A  A 

V z  2 ;  —  .  ^X  X  \  Vu  11 

A  +  2 fi  ,x  A  +  2/r  y'y 

Vy,z  = 

Vx,z  =  —Vz,x- 


(26) 

(27) 

(28) 


Taking  these  equations  into  the  stress-strain  relation,  the  real  equations  to  update  axx  and  oyy  at  the 
free  surface  are 


&xx,t  —  (A  T  2fi^vx^x  T  A Vy^y  -f-  A 
®yy,t  —  Xvx,x  T  (A  +  2/i)r+;^  T  A 


A 

77  V  nr  x 

A  +  2  fi  ’ 

A 

V  02  x 

A  +  2/r  x’x 


A  +  2  fj,Vy'y_ 


\  +  2nVy’y 


(29) 

(30) 


To  effectively  absorb  the  surface  wave,  the  ADE  CFS-PML  equations  of  axx  and  ayy  at  the  inter¬ 
section  of  the  free  surface  and  the  PML  should  be  derived  from  the  above  equations  as 

_  _  (A  +  2+)2  -  A2  1  ,  2A ijl  1  (A  +  2/r)2  -  A2  1  TXX  2A+  1  TXV 

xx't_  A  +  2/r  pxVx’x+  \  +  2^(3yVy’y  X  +  2fi  A  +  2^/^ 


(31) 

2A/i  1  ,  (A  +  2fi)2  -  A2  1  2A/u  1  TXX  (A  +  2^f  -  A2  1  TXJ/ 

<7w,t  \  +  2^pxVx’x+  A  +  2M  PyVy’y  \  +  2vPxx  A  +  2M  PyV' 

(32) 


The  above  equations  in  the  correction  form  to  the  PML  equations  is  shown  in  Appendix  A. 


In  the  FDTD  method,  there  are  two  possible  locations  for  the  free  surface,  coincident  with 
the  normal  stress  components  or  with  the  shear  stress  component  axz  (Gottschammer  and  Olsen, 
2001;  Kristek  et  ah,  2002).  As  the  numerical  tests  presented  in  the  following  section  show,  the 
above  modified  ADE  CFS-PME  equations  are  crucial  to  effectively  absorb  surface  waves  if  the  free 
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surface  coincides  with  the  normal  stress  components.  If  the  free  surface  coincides  with  axz,  the 
free  surface  boundary  condition  is  implemented  as  in  the  physical  domain  by  setting  axz  and  ayz 


to  zero  at  the  free  surface  and  no  special  treatment  is  needed  in  the  PML  layer.  In  the  nonstaggered 


finite-difference  method,  our  numerical  tests  show  that  the  modified  PML  equations  3 1  and  32  are 
essential  for  the  stability  of  the  solution  if  the  velocity  component  free  surface  boundary  condition 
(equation  26)  is  applied. 

THE  ROLE  OF  THE  SCALING  FACTOR  fi 

In  the  context  of  electromagnetic  wave  modeling,  Petropoulos  (2000)  stated  that  f3x  “allows  for 
bending  the  waves  towards  the  normal  direction  into  the  layer  by  increasing  the  corresponding 
component  of  the  wave  vector  while  leaving  the  tangential  component  unchanged;  at  the  same 
time,  it  increases  the  damping  in  the  normal  direction  into  the  layer”  with  a  stretching  function 
sx  =  /?x(l  +  n  d\i.)-  The  role  of  parameter  fix  has  not  been  thoroughly  understood  in  seismic  wave 
modeling.  For  example,  a  constant  (3  =  1  was  used  in  Festa  and  Vilotte  (2005)  and  Festa  et  al. 
(2005).  fi  did  not  significantly  affect  the  results  in  Komatitsch  and  Martin  (2007)’s  numerical  tests, 
but  played  an  important  role  in  Drossaert  and  Giannopoulos  (2007a,b)’s  numerical  tests. 

To  understand  the  role  of  fix,  we  rewrite  the  CFS  stretching  function  (equation  12)  as 


(33) 


The  complex  stretching  of  the  velocity-stress  equation  1-2  with  the  stretching  function  equation  33 
is  equivalent  to  using  a  stretching  function 


dx/  fix 


(34) 
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to  stretch  the  following  equations 


dv  1  /  1  „  „ 

"777  —  ~  (  ~R~ ^1 

dt  p  \fix 


0, 


da 


=  c  : 


1 


“  Vo;  +  V II  )  V  + 


Vz  +  V II  V 


-i  T ' 


(35) 


dt  ‘  [  \fix 

The  solutions  of  equation  35  in  a  homogenous  isotropic  medium  is  the  same  as  the  solution  of  the 

velocity-stress  equations  1-2  in  a  transversely  isotropic  medium  ( x  is  the  symmetry  axis)  that  has 

the  stiffness  tensor 

0  0  0 

0  0  0 


c  = 


A+2/i 

~w 


\_ 

Px 


A 

px 


\ 


fix 


V 


A 

Px 

0 

0 

0 


A  T  2/r  A 

A  A  +  2p  0  0  0 

0 


0 

0 

0 


p  0  0 

0 


m 


0  0 
0  0  0 


t1 

M) 


(36) 


Thus  (3X  introduces  an  anisotropy  that  reduces  the  phase  velocity  of  the  waves  propagating  normal 


to  the  PML  layer  to  The  plane  wave  solution  of  equation  35  can  be  written  as, 


P-. 


u  =u0eia;(t-sx), 


(37) 


where  uo  is  the  polarization  vector,  s  =  is  the  slowness  vector  which  is  normal  to  the  wavefront 
surface,  C  is  the  phase  velocity  depending  on  the  propagation  direction.  For  simplicity,  we  considier 
waves  propagating  in  (x,  z)-plane.  Thus  equation  37  becomes 


u  =uoe 


(38) 


where  6  is  the  angle  between  the  normal  direction  of  the  wavefront  and  x-axis.  Taking  equation  34 
into  equation  8,  we  can  get  the  stretched  x 


x  =x  + 


f 

Jo 


Otxdx/  fix 

ofi.  +  L 02 


d  r)  —  i 


ojdx/ fix 

+  u;5 


d?y. 


(39) 
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Taking  above  equation  into  equation  38,  we  get  the  plane  wave  solution  of  the  CFS-PML 


u  =  Une  V  c  c  )e  C  JO 


x  Q->dx  / 0x °^x 


ax+u 


dr] 


cos  6 

c 


rx 

Jo 


ydx/ [3X  dr I 


(40) 


D 

For  waves  propagate  normal  to  the  PML  layer  (phase  velocity  C  =  ^  ,  and  cos  9  =  1,  assuming 
a  constant  j3x ),  the  damping  term  D  becomes  exp(— ^[Jq  Jf+ji dx  dr/]),  in  which  j3x  is  canceled 
out.  This  indicates  that  for  normal  incident  waves,  j3x  only  decreases  the  wave  propagation  speed 
without  increasing  damping.  This  is  still  true  for  a  varying  fJx  profile,  since  each  single  layer 
within  the  PML  can  be  taken  as  having  a  constant  /3.  We  validate  this  by  calculating  the  waveform 
in  the  PML  with  different  (3X  values  for  a  normal  incident  plane  wave  using  the  FDTD  method. 
We  set  ax  =  0  and  use  polynomial  scaled  /3X  and  dx.  Figure  2a  shows  the  maximum  amplitude 
for  different  j3  values  (3o  =  1  and  3o  =  9)  at  each  grid  point  inside  a  25-cell-thick  PML  layer. 
The  maximum  amplitude  distributions  are  nearly  identical  for  different  /?,  verifying  that  /3  does 
not  increase  damping  directly.  Figure  2b  is  the  comparison  of  velocity  at  the  grid  point  15 -cell 
into  the  PML,  in  which  the  wave  using  /So  =  9  arrives  later  than  that  using  /So  =  1  due  to  the 
reduced  wave  speed  by  (3X.  We  note  that  the  stretching  function  in  Petropoulos  (2000)  is  of  a  form 
sx  =  3X(1  +  a  d+ibJ),  whose  damping  term  contains  the  product  of  3X  and  dx  in  exponent.  Therefore 
fjx  in  Petropoulos  (2000)  affects  damping. 


For  grazing  incident  waves,  cos  9  is  very  small  without  using  3X.  Thus  the  damping  term  D  in 
equation  38  may  not  be  small  enough  to  efficiently  damp  the  wave  propagating  sub-parallel  to  the 
PML  layer  and  causes  strong  spurious  reflections.  Because  of  the  anisotropy  generated  by  {3X,  the 
normal  direction  of  the  wavefront  is  bended  towards  the  normal  direction  of  the  PML  layer,  thus 
cos  9  in  the  damping  term  is  increased  and  the  absorption  is  improved.  It  should  be  noted  that  3x 
can  not  be  canceled  out  in  the  damping  term  for  oblique  incident  waves  because  the  phase  velocity 
is  a  function  depending  on  several  components  of  the  stiffness  tensor. 
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NUMERICAL  TESTS 


To  demonstrate  the  ADE  CFS-PML,  we  present  results  of  three  numerical  tests.  One  is  performed 
in  a  thin  homogeneous  3D  slab  surrounded  by  PMLs  on  all  surfaces  to  investigate  the  optimal 
values  of  the  CFS  parameters.  In  the  second  test,  we  add  a  free  surface  to  the  thin  slab  model  to 
validate  the  modified  PMF  equations  that  take  into  accout  the  free  surface  boundary  condition  in 
absoiption  of  surface  waves.  We  use  the  FDTD  method  in  the  first  two  tests.  In  the  third  test,  we  use 
the  nonstaggered  finite-difference  method  with  the  4th-order  Runge-Kutta  scheme  in  a  two-layered 
model  with  a  3D  topography  to  validate  the  ADE  CFS-PMF  implementation  in  a  higher-order  time 
marching  scheme. 

Before  showing  the  numerical  test  results,  let  us  briefly  review  the  choice  of  the  parameter 
values  in  the  CFS  function.  According  to  Roden  and  Gedney  (2000),  d,  (3  and  a  should  be  spatially 
scaled  in  the  PMF  layer  to  reduce  reflection  errors  in  the  discrete  space,  d  usually  is  zero  at  the 
PMF-interior  interface  and  maximum  at  the  exterior  boundary.  f3  is  equal  to  1  at  the  PMF-interior 
interface  and  maximum  at  the  exterior  boundary,  a  is  scaled  in  a  reversed  manner  that  the  value  is 
maximum  at  the  PMF-interior  interface  and  gradually  reduced  to  zero  at  the  exterior  boundary  to 
absorb  low  frequency  waves. 

In  this  paper,  we  choose  the  commonly  used  p-order  polynomial  scaling  function  (Gedney, 
1998;  Roden  and  Gedney,  2000), 


(41) 


(42) 


(43) 


where  do  and  (3q  are  the  maximum  value  of  d  and  3  at  the  exterior  boundary  (the  optimal  (3q  is 
discussed  in  the  next  section)  ,  op  is  the  maximum  value  of  a  at  the  PMF-interior  interface.  L  is 
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the  width  of  the  PML  layer,  x  is  the  distance  to  the  PML-interior  interface.  p,j,  pp  and  pa  typically 
range  from  1  to  4,  and  2  is  commonly  used.  We  use  pd  =  2  and  pp  =  2.  We  choose  pa  =  1  since 
Martin  et  al.  (2008b)  showed  the  linear  variation  of  a  gets  a  more  pronounced  decay  of  energy. 


The  maximum  value  of  d  is  obtained  through  a  theoretical  reflection  coefficient  relation  (Collino 
and  Tsogka,  2001), 


do  — 


(; Pd  +  1  )cp 
2  L 


In  R, 


(44) 


where  cp  is  the  compressional  wave  speed,  and  R  is  the  theoretical  reflection  coefficient  for  a  normal 
incident  plane  P  wave  with  a  Dirichlet  condition  (v  =  0  and  a  =  0)  at  the  exterior  boundary  of 
the  PML  layer.  The  trade-off  between  a  small  and  a  large  do  is  that  a  too  small  do  can  not  damp 
the  wavefield  enough  in  the  PML  layer  and  strong  spurious  reflections  from  the  Dirichlet  boundary 
will  propagate  back  to  the  physical  domain,  and  a  too  large  do  leads  to  spurious  reflections  caused 
by  the  dispersion  (Collino  and  Tsogka,  2001).  Thus  the  optimal  do  depends  on  the  PML  thickness 
L  and  the  theoretical  reflection  coefficient  II.  Collino  and  Tsogka  (2001)  proposed  R  for  different 
width  of  the  PML  layer  as  R5  =  0.01,  Rw  =  0.001  and  R 20  =  0.0001,  where  the  superscript 
denotes  the  width  of  the  PML  region  in  the  amount  of  grid  spacing.  Marcinkovich  and  Olsen 
(2003)  approximated  the  above  relation  by  a  third  order  polynomial  function  to  get  the  R  value  for 
an  arbitrary  PML  width  between  1  and  20  layers.  We  approximate  the  relation  by  a  logarithmic 
equation, 


(  R  \  _  log10(R20/Rw)  1  /AC 
°gl0U10J  logi0(20/10)  10610  U, 


(45) 


Thus  R  for  an  arbitrary  PML  layer  width  of  N  cell  size  can  be  expressed  as 


,  „  /m  _  log10(JV)  -  1  0 

l0gl»(-R)  -  log10(2)  “ 3' 


(46) 


Equation  44  is  optimized  for  normal  incident  waves.  For  oblique  incident  waves,  a  larger  do  than 
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that  from  equation  44  is  needed  to  obtain  optimal  damping.  From  numerical  tests,  we  find  that 
increasing  do  to  2-3  times  the  value  from  equation  44  yields  better  solutions  for  grazing  incidence. 

Festa  and  Vilotte  (2005)  recommended  op  to  be  nfc,  where  fc  is  the  dominant  frequency  of  the 
source  time  function.  In  our  numerical  tests,  different  values  of  op  (1.07r/c,  1.57 r/c  and  2.0irfc)  do 
not  make  a  significant  difference.  Thus  op  =  "  fc  seems  a  good  choice  in  practice. 

Thin  Homogeneous  3D  Slab  Surrounded  by  PMLs 

This  test  compares  the  performance  of  the  ADE  CFS-PML  with  the  standard  PML  (by  setting 
(3o  =  1  and  op  =  0)  for  body  waves  impinging  the  PML  interface  at  a  near-grazing  angle  and 
shows  how  the  values  of  do,  /3q  and  o:q  affect  the  absorption. 

Figure  3  shows  the  computational  configurations.  The  source  time  function  is  a  Ricker  wavelet 
with  a  center  frequency  of  5  Hz  and  a  time  delay  of  0.6  second.  The  time  step  is  1.5  ms  and  the 
total  time  step  is  5000,  resulting  in  a  time  window  of  7.5  s,  in  which  the  reflected  S  wave  can  travel 
through  the  computational  domain.  The  value  of  do  is  obtained  through  equation  44.  Different 
f3o  and  different  op  (0, 0.57 r/c,  1.07r/c,  1.57t/c  and  2.07 r/c)  are  tested.  The  result  of  a  large  model 
(double  the  size  along  x  and  y  and  twelve  times  larger  along  z  with  a  20-cell-thick  PML)  is  used  as 
the  reference. 

Ligure  4  shows  the  snapshots  of  the  vz  component  at  3.75  s.  P  wave  has  propagated  outside 
the  computational  domain  and  the  dominant  phase  is  S  wave.  The  standard  PML  (Ligure  4b)  shows 
stronger  spurious  evanescent  waves  at  the  positive  z  face,  o-q  =  7r/c  (Ligure  4c)  can  eliminate  such 
kind  of  spurious  waves,  but  a  spurious  reflected  wave  immediately  after  the  S  phase  can  be  clearly 
seen.  Increasing  oo  to  27 r/c  does  not  significantly  affect  the  result  (not  shown).  Using  3q  =  10 
in  addition  to  ao  =  7r/c  (Ligure  4d)  reduces  the  spurious  reflected  wave  and  yields  better  solution 
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for  grazing  incidence.  Increasing  do  to  three  time  larger  (Figure  4e)  can  dramatically  improve 
the  absorption  of  grazing  incident  waves,  but  causes  a  stronger  spurious  reflection  of  waves  at  near¬ 
normal  incidence  at  the  positvie  y  face.  The  above  results  indicate  that  the  choices  of  the  parameters 
Px,  ax  and  do  are  important  for  absorption  of  grazing  incident  waves.  A  thicker  PML  (20-cell-thick, 
Figure  4f)  can  further  absorb  the  normal  incident  waves. 

Figure  5  shows  the  comparison  of  the  vz  time  series  between  the  PML  models  and  reference 
model  at  the  receiver  indicated  in  Figure  3.  The  amplitude  of  S  wave  in  both  the  standard  PML 
(Figure  5a)  and  the  ADE  CFS-PML  with  o-o  =  7 r/c  and  Pq  =  1  (Figure  5c)  is  larger  than  that 
of  the  reference  model.  This  larger  amplitude  S  wave  consists  of  the  incident  S  wave  and  the 
spurious  reflected  S  wave,  which  propagates  in  the  nearly  same  direction  as  the  incident  S  wave. 
Figure  5c  suggests  that  though  a  can  eliminate  the  spurious  evanescent  waves,  it  does  not  increase 
the  damping  of  waves  at  grazing  incidence,  thus  the  spurious  reflected  body  waves  are  still  strong. 
This  is  reasonable  from  the  damping  term  expression  (without  fix)  exp(— f()r  dx  d//)  that 
a  does  not  significantly  affect  the  incident  angle  0,  thus  cos  9  is  very  small  for  near-grazing  incident 
waves  no  matter  what  value  of  a  is.  On  the  contrary,  0X  bends  the  wavefront  towards  the  normal 
direction  of  the  PML,  thus  decreases  the  slowness  angle  9  and  increases  the  damping  of  waves  at 
near-grazing  incidence.  This  effect  can  be  seen  from  a  snapshot  comparison  of  0  =  1  (Figure  6a) 
and  p  =  10  (Figure  6b).  This  is  the  reason  that  Figure  5e  (cto  =  tt/c  and  Pq  =  10)  shows  the 
amplitude  of  S  wave  is  close  to  that  of  the  reference.  But  the  increased  damping  due  to  bending 
is  not  large  enough  for  P  wave,  so  a  discrepancy  in  P  wave  can  be  seen  in  Figure  5f.  A  three 
times  larger  do  (with  0o  =  7,  ao  =  vr  /',,)  significantly  improves  the  absorption  of  waves  at  grazing 
incidence  and  has  nearly  perfect  fits  for  both  P  and  S  waves  (Figure  5g-5h). 

To  quantitatively  measure  the  error,  we  use  the  following  integral  as  a  local  error  at  the  position 
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x. 


£l(x)  = 


T,t  bpml(x^)  -  ^ref(x,t)| 


(47) 


£t  l^ref  (Xj  f  )  | 

where  fpmi(x,  t)  is  the  simulated  wavefield  in  the  PML  model  and  nref(x,  t)  is  the  counterpart  in  the 
reference  model.  The  local  error  distribution  (Figure  7)  confirms  the  conclusions  from  the  waveform 
comparisons,  that  a  (Figure  7b,  ao  =  vr/c,  Pq  =  1)  eliminates  the  spurious  evanescent  waves  in 
the  standard  PML  (Figure  7a)  and  reduces  the  local  error  in  the  upper  part  of  the  computational 
domain,  but  the  maximum  error  is  still  on  the  same  order  because  the  amplitude  of  the  reflected  S 
wave  is  similar;  p  (Figure  7c,  cro  =  vr/c,  Pq  =  10)  improves  the  absorption  of  the  grazing  incident 
waves  and  the  local  error  is  much  smaller;  a  three  time  larger  do  efficiently  increases  damping  of 
the  grazing  incident  waves  and  reduces  the  local  error  (Figure  7d,  ao  =  vr  fc,  Pq  =  7). 


We  use  an  averaged  local  error  in  a  subgrid  located  in  the  xz  plane  through  y  =  3000m  as  an 

indicator  of  the  global  error  of  the  PML  with  a  particular  Pq  and  ao, 

1  580  301  61 

eg  =  ^  E  E  (48) 

V  i=20  j=301  fc=20 

where  N  is  the  total  number  of  the  grid  points  of  the  subgrid,  i,  j.  k  are  the  grid  index  in  the 
subgrid.  The  subgrid  is  10-cell  away  from  the  PML-interior  interfaces.  Figure  8a  shows  the  contour 
plot  of  the  global  error  as  a  function  of  Pq  and  ao-  It  can  be  seen  that  ao/(7r/c)  between  1  to 
2  does  not  affect  the  overall  accuracy.  The  optimal  value  of  Pq  is  around  10  ~  12  in  this  case. 
To  investigate  the  optimal  Pq  variation  with  respect  to  the  PML  thickness  N  and  the  grid  spacing 
(or  wave  speed,  or  source  frequency),  we  perform  a  series  numerical  calculations  with  different  N 
and  cs.  The  global  error  with  respect  to  Pq  in  Figure  8b  shows  that  the  optimal  Pq  is  insensitive 
to  the  PML  thickness  N,  but  inreases  when  cs  increases  (equivalent  to  grid  spacing  decreases,  or 
source  frequency  decreases).  Theoretically,  a  larger  P  can  bend  waves  more  efficiently,  but  a  very 
large  P  can  make  the  wavelength  too  short  to  be  resolved  by  the  numerical  scheme,  which  can 
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cause  strange  numerical  spurious  reflections.  Figure  8c  and  8d  show  the  normalized  points  per 
dominant  wavelength  (=  PPW/PPWo,  where  PPW  is  points  per  dominant  wavelength,  PPWo  is  the 
minimal  PPW  requirement  of  the  numerical  scheme.  For  the  4th-order  staggered  finite-difference, 
PPWo  =  6)  distribution  inside  the  PML  layer  at  the  center  frequency  for  the  optimal  Pq  in  Figure  8b 
(/3o  =  11  for  c.5  =  1600  nr/s,  and  (3q  =  15  for  cs  =  2300  nr/s).  The  normalized  points  per  dominant 
wavelength  at  the  exterior  boundary  is  around  0.5  in  difference  cases.  Thus  we  can  estimate  the 
optimal  fio  through 


Po 


C 

0.5PPW0  A hfc' 


(49) 


It  should  be  noted  that  we  use  the  same  parameter  values  in  all  PMLs  on  the  six  faces  in  this 
numerical  test.  In  practice,  one  can  use  larger  /3o  and  do  only  for  a  particular  PML  layer,  for  which 
grazing  incidence  is  important. 


Thin  3D  Slab  With  the  Free  Surface 

We  change  one  face  of  the  thin  slab  in  the  previous  test  to  the  free  surface  to  demonstrate  the 
absorption  of  surface  waves  of  the  ADE  CFS-PML  and  the  effect  of  the  modified  PML  equations 
that  take  into  account  the  free  surface  boundary  condition. 

Figure  9  shows  the  computational  configurations.  The  vertical  force  is  located  at  (200,  3000,  795) 
nr,  half  grid  below  the  free  surface.  All  other  parameters  are  the  same  as  the  model  in  Figure  3.  The 
free  surface  boundary  condition  is  implemented  using  AFDA  technique  (Kristek  et  ah,  2002),  in 
which  a  compact  finite-difference  operator  and  biased  finite-difference  operators  are  used  to  calcu¬ 
late  derivatives  with  respect  to  z  axis  at  grid  points  near  the  free  surface. 

Figure  10  is  the  snapshot  of  the  vz  component  at  3.75  s.  The  free  surface  coincides  with  the 
location  of  the  normal  stress  components.  Because  waves  do  not  impinge  the  PML  at  near-grazing 
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incidence,  the  standard  PML  (Figure  10b)  and  the  ADE  CFS-PML  (Figure  10c,  (3q  =  1  and  op  = 
7 r/c)  both  efficiently  absorb  surface  waves.  For  the  same  reason,  increasing  (3q  (not  shown  here) 
does  not  make  significant  difference.  If  the  modified  ADE  CFS-PMF  equations  31  and  32  which 
incorporate  the  free  surface  boundary  are  not  used,  a  relatively  strong  reflected  surface  phase  can 
be  clearly  seen  (Figure  lOd).  If  the  free  surface  is  located  at  the  axz  component,  the  amplitude  of 
the  reflected  wave  is  smaller  than  Figure  lOd,  but  stronger  than  Figure  10b.  Thus  the  absorption  of 
surface  waves  when  the  free  surface  coincident  with  the  axz  component  is  not  as  efficient  as  the  free 
surface  coincident  with  the  normal  stress  components. 

Two- layered  Medium  With  Free  Surface  Topography  in  a  Narrow  Vertical  Slice 

In  this  test,  we  validate  the  ADE  CFS-PMF  implementation  in  the  4th-order  Runge-Kutta  scheme, 
and  demonstrate  the  optimal  performance  of  the  ADE  CFS-PMF  in  a  narrow  vertical  slice  to  simu¬ 
late  wave  propagation  in  a  common  situation  involving  seismic  profiles. 

The  medium  has  two  layers  (Figure  11)  separated  at  z  =  —8km.  Boundary-conforming  grid  is 
used  to  conform  grid  with  the  surface  topography.  The  horizontal  grid  spacing  is  50m  and  a  varying 
grid  spacing  (33m  ~  60m)  is  used  vertically.  The  source  time  function  is  a  Ricker  wavelet  with  a 
center  frequency  of  1.5  Hz  and  a  time  delay  of  1  s.  The  value  of  do  is  calculated  from  equation  44, 
where  the  cp  of  the  first  layer  (3000  m/s)  is  used  for  the  PMF  layers  at  x  and  y  faces,  and  cp  of 
the  lower  layer  (5000  m/s)  is  used  for  the  bottom  PMF  layer.  The  value  of  op  is  7 r/c.  The  optimal 
/?o  is  estimated  to  be  8  from  equation  49  using  the  surface  wave  speed  of  the  first  layer.  We  also 
simulate  with  (30  =  1  for  comparison  with  the  optimal  (3q  value.  Numerical  solution  of  a  larger 
model  (double  the  length  and  depth,  five  times  the  width,  a  25 -cell-thick  PMF)  is  used  as  reference. 

The  difference  between  the  ADE  CFS-PMF  (Figure  12b,  ao  =  7r/c  and  ;3q  =  1)  solution  and 
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that  of  the  reference  (Figure  12a)  is  so  small  that  it  is  hard  to  observe  from  the  snapshots  when 
the  color  is  scaled  to  1%  of  the  maximum  amplitude.  Thus  we  compare  the  vz  time  series  at  the 
receiver  (triangle  in  Figure  11).  Figure  13a  shows  the  comparison  between  the  solution  of  the 
reference  (solid  line),  PML  with  [3q  =  1  (dash-dotted  line)  and  PML  with  /5o  =  8  (dotted  line), 
ao  =  vr/c  is  used  for  both  PML  solutions.  The  fit  is  nearly  perfect  at  this  scale.  Figure  13b  shows 
the  comparison  in  the  time  window  13  ~  21  s  with  amplitude  scaled  to  0.2%  of  that  in  Figure  13a. 
We  note  that  even  at  this  scale,  the  solution  of  /3q  =  8  (dotted  line)  overlaps  with  the  reference, 
suggesting  that  the  ADE  CFS-PML  with  the  optimal  (3  matches  the  reference  for  this  complex 
waveform  containing  scattered  waves  from  the  topography  and  multiple  reflections  from  the  inner 
interface.  To  further  highlight  the  difference,  Figure  13c  shows  the  vz  difference  between  the  PML 
with  fio  =  1  and  the  reference  (nf0_1  —  vrze^)  and  that  between  the  PML  with  =  8  and  the 
reference  (vz~S  —  vrz^).  The  y- axis  scale  is  around  0.5%  of  that  in  Figure  13a.  It  is  clear  the 
ADE  CFS-PML  can  efficiently  absorb  waves  in  this  narrow  two-layered  topographic  model.  The 
spurious  waves  do  not  significantly  affect  the  reflections  from  the  inner  interface.  The  optimal  /3q 
yields  a  better  solution  with  the  amplitude  of  the  spurious  waves  smaller  than  that  of  (3q  =  1,  at  a 
level  below  0.1%  of  the  maximum  amplitude.  This  numerical  test  indicates  that  we  can  use  a  narrow 
slice  mesh  with  the  ADE  CFS-PML  to  simulate  full  wave  propagation  in  complex  structure  models 
in  a  common  situation  involving  seismic  profiles. 

CONCLUSION 

We  introduce  an  efficient  unsplit-field  complex  frequency-shifted  perfectly  matched  layer  imple¬ 
mentation  using  auxiliary  differential  equations  in  seismic  wave  modeling.  We  derive  the  complete 
ADE  CFS-PML  equations  for  the  velocity-stress  equations.  Because  the  ADE  CFS-PML  equations 
are  first-order  partial  differential  equations,  it  can  be  solved  by  the  some  numerical  scheme  used 
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in  the  inner  domain  no  matter  what  time  marching  scheme  is  used.  We  demonstrate  this  by  im¬ 
plementing  the  ADE  CFS-PML  in  the  2nd-order  leap-frog  scheme  and  the  4th-order  Runge-Kutta 
scheme. 

We  analyze  the  role  of  the  scaling  factor  (3  in  the  CFS-PML  and  demonstrate  that  it  introduces 
an  intermediate  transformation  of  the  velocity-stress  equations  into  an  anisotropic  medium  that  has 
a  lower  wave  speed  normal  to  the  PML  layer.  For  normal  incident  waves,  it  only  delays  the  waves 
propagating  into  the  PML  layer  without  increased  damping.  For  oblique  incident  waves,  it  can  bend 
the  waves  towards  the  normal  direction,  thus  increases  cos  9  in  the  damping  term  and  improves 
the  absorption.  Numerical  tests  indicate  that  the  frequency  shifted  factor  a  only  eliminates  the 
evanescent  waves  caused  by  near-grazing  incident  waves,  but  does  not  significantly  improve  the 
absoiption  of  the  reflected  body  waves.  (3  is  more  efficient  and  important  in  absorbing  grazing 
incident  waves. 

Our  numerical  tests  indicate  that  the  incompatibility  between  the  PML  equations  and  the  free 
surface  boundary  condition  implementation  can  cause  an  unstable  wave  originating  from  the  free 
surface  in  the  PML  layer  in  the  nonstaggered-grid  finite-difference.  Thus  we  derive  the  PML  equa¬ 
tions  taking  the  free  surface  boundary  condition  into  account.  Numerical  tests  validate  that  these 
modified  PML  equations  solve  the  instability  problem  in  the  non-staggered  finite-difference  method, 
and  improve  the  efficiency  of  absorbing  surface  waves  in  the  staggered  finite-difference  method 
when  the  free  surface  coincides  with  the  normal  stress  components. 

We  perform  a  series  of  numerical  tests  to  validate  the  ADE  CFS-PML  implementation  both  in 
the  FDTD  method  and  the  nonstaggered-grid  finite-difference  method  using  the  4th-order  Runge- 
Kutta  scheme.  We  investigate  the  different  values  of  the  CFS  parameters  on  the  efficiency  of  absorp¬ 
tion  based  on  quantitative  error  criteria.  Our  numerical  tests  indicate  that  the  optimal  (3  reduces  the 
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points  per  dominant  wavelength  at  the  exterior  boundary  to  three  for  the  4th-order  finite-difference 
operator,  around  half  the  value  required  by  the  numerical  scheme.  Different  values  of  ao  (1 .07 r/c, 
1.57t/c  and  2.0irfc)  do  not  make  a  significant  difference.  Thus  «o  =  7r/c  seems  a  good  choice  in 
practice.  For  grazing  incidence,  adjusting  only  parameters  a  and  /3  is  not  enough,  one  also  needs  to 
increase  the  maximum  value  do  of  the  damping  profile  d  to  get  optimal  absoiption. 

Numerical  test  of  a  two-layered  model  with  a  surface  topography  in  a  narrow  vertical  slice 
indicates  that  the  ADE  CFS-PML  yields  accurate  results  compared  with  the  reference  in  a  much 
larger  and  thicker  model,  thus  making  it  possible  to  carry  out  efficient  3D  forward  simulation  in  a 
common  situation  involving  seismic  profiles  in  seismic  exploration. 
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APPENDIX  A 

ADE  CFS-PML  EQUATIONS  FOR  THE  ELASTIC  WAVE  EQUATIONS 

ADE  CFS-PML  equations  for  the  velocity-stress  formulation  of  the  elastic  wave  equations  are 


&XX,t  (A  +  2 /l)VXix  A Vy^y  +  A vZjZ 

+  (A  +  2/i) 


'  1 

r  i 

'  1 

Jx  -  1 

Vx,x  +  A 

i 

i — i 

_ 1 

vv,y  +  ^ 

Jz~l. 

-  (A  +  2»)±-VZ  -  X^VV  - 
Dy 


"Py  "  "Pz 


&yy,t  — XVx.x.  A  (A  A  ‘2^-ljVy^y  A  XvZ)Z 


AA 


Px 


Vx,x  A  (A  A  2/i) 


'  1 

'  1 

.Py 

-  1 

vy,y  +  ^ 

Jz~  . 

&zz,t  —XvX}X  A  A vVtV  A  (A  A  2 /i)vZtZ 


AA 

'  1 

- 1 

^x,x  H-  A 

2-  -i 

.Px 

.Py 

vv,v  +  (2v  A  2/i) 


Pz 


- 1 


Vyz,t  =P>  ( Vy,z  +  VZtV)  +  /i 
&xz,t  =Pj  ( Vx,z  H"  Vz, rr)  H-  P 
&xy,t  =P  (Vx,y  ^ y,x )  “1“  P 
pVx,t  — &xx,x  “1“  &xy,y  H-  &xz, 


Px 


P^y,t  — &xy,x  “1“  &yy,y  H-  &yz,z 


k~\ 

h~\ 


wy,z 


Py  1. 

.Px 


-  1 


1 

Px 


-  1 


Py 


®XV,V 


T,~\ 


'  1 

'  l 

'  1 

A 

.Px 

- 1 

&xy,x  H"~ 

.Py 

-  1 

ayy,y  + 

Jz  ~  . 

pvz,t  ~&xz,x  +  0^2,  y  +  ^22,2 


'  1 

r  i 

'  1 

A 

Jx~  . 

®XZ,X  “t“ 

i 

i— 1 

ayz,y  + 

Jz  ~  . 

rr _ Tx _ Ty _ Tz 

uxz,z  q  -1  xx  a  -1  xy  o  ±xz'> 

Px  Py  Pz 


T _ L  Tx _ — Tv _ —  Tz 

Vziz  p  xy  p  zyy  p  -Lyzt 


cr~~  ~ - Tx - Tv - Tz 

Z  ’Z  Px  XZ  Py  VZ  Pz  Z2’ 


(A- 1 ) 


(A-2) 


(A- 3) 
(A-4) 
(A-5) 
(A-6) 

(A-7) 

(A-8) 

(A-9) 


26 


A-2  7 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 


GEOPHYSICS 


with  the  auxiliary  differential  equations  for  the  memory  variables  damping  along  x. 
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At  the  free  surface  in  the  PML  layer,  the  PML  equations  of  <Txx,t  and  (Jyy.t  taking  the  free  surface 
boundary  condition  into  account  are 
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LIST  OF  FIGURES 

1  A  3D  computational  domain  surrounded  by  PML  layers  on  each  face.  Dashed  lines  mark 
the  PML-interior  interfaces.  Along  the  edge  of  the  computational  domain,  two  directional  PMLs 
exist.  At  the  corner,  three  directional  PMLs  co-exist. 

2  (a)  Maximum  amplitude  distribution  for  different  P  values  (Pq  =  1,  solid  line  and  Po  =  9, 
dotted  line)  for  a  25 -cell-thick  PML  of  a  normal  incident  plane  wave,  a  =  0  and  an  optimal  cL 
are  used,  (b)  Velocity  seismogram  at  the  grid  point  15-cell  into  the  PML  for  P o=l  (solid  line)  and 
/?o  =  9  (dotted  line).  The  wave  using  Po  =  9  arrives  later  than  that  using  Po  =  1  due  to  the  reduced 
wave  speed  by  /3X. 

3  A  thin  homogeneous  3D  slab  surrounded  by  10-cell-thick  PMLs  on  all  six  faces  (three 
cross-sections  are  shown).  The  PML-interior  interfaces  are  indicated  by  dashed  lines.  The  medium 
has  a  compressional  wave  speed  cp  =  3200  m/s,  shear  wave  speed  cs  =  1600  m/s,  and  density 
p  =  2200  kg/m3.  The  grid  spacing  is  10  m,  the  total  grid  number  is  600  x  600  x  81,  and  the  size 
of  the  model  is  5990  x  5990  x  800  m.  The  star  symbol  shows  the  location  of  the  vertical  force 
(200,  3000,  600)  m,  ten  cells  away  both  from  the  PML-interior  interface  at  the  negative  x  face  and 
from  that  at  the  positive  2;  face.  The  triangle  symbol  indicates  the  receiver  (5790,  3000,  600)  m,  ten 
cells  away  from  both  the  PML-interior  interface  at  the  positive  x  face  and  that  at  the  positive  2  face. 

4  Snapshots  of  the  vz  component  of  the  thin  slab  model  (Figure  3)  at  3.75  s.  The  color  is 
scaled  to  0.1%  of  the  maximum  absolute  amplitude.  White  lines  are  the  PML-interior  interfaces, 
(a)  Reference  solution,  (b)  Standard  PML.  (c)  oq  =  vr/c  and  Pq  =  1.  (d)  oq  =  vr fc  and  Pq  =  10. 
(e)  O'o  =  7r/c,  Pq  =  7  and  a  three  times  larger  do-  (f)  «o  =  7r/c  and  Pq  =  11  for  a  20-cell-thick 
PML. 

5  The  vz  seismograms  at  the  receiver  shown  in  Figure  3.  The  right  column  shows  the  same 
seismograms  except  that  the  y  axis  is  scaled  to  1%  of  that  of  the  left  column  to  highlight  the  dif- 
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ference.  The  first  arrival  in  the  right  column  is  P  wave,  which  is  nearly  invisible  in  the  left  column. 

(a) -(b)  Standard  PML.  (c)-(d)  To  =  1  and  a0  =  tt/c-  (e)-(f)  To  =  10  and  a0  =  7 r/c.  (g)-(h)  To  =  7, 
o:o  =  7 r/c  and  a  three  times  larger  do- 

6  Snapshots  of  vz  on  the  cross-section  in  the  middle  of  y  axis  at  3.75  s.  The  color  is  scaled  to 
the  same  maximum  absolute  amplitude.  Green  lines  are  the  PML-interior  interfaces,  (a)  «o  =  7r/c 
and  /3 o  =  1-  The  amplitude  of  S  wave  in  the  upper  PML  layer  is  even  larger  than  that  in  the  physi¬ 
cal  region,  which  indicates  a  does  not  efficiently  absorb  this  near-grazing  propagating  body  wave. 
The  distribution  of  the  peak  amplitude  is  wider  than  Figure  6b,  which  indicates  the  reflected  wave 
propagates  into  the  physical  domain,  (b)  op  =  7r/c  and  To  =  10.  The  wave  front  is  bended  towards 
the  normal  direction  inside  both  the  upper  PML  layer  and  the  bottom  PML  layer.  The  amplitude  of 
S  wave  exhibits  damping  in  the  PML  layers. 

7  Local  error  distribution.  The  color  is  scaled  to  0.5.  White  lines  show  the  PML-interior 
interfaces,  (a)  Standard  PML.  (b)  To  =  1  and  a?o  =  ^ t/c-  (c)  To  =  10  and  op  =  71  fr-  (d)  To  =  7, 
oo  =  7 r/c  and  a  three  times  larger  do- 

8  (a)  Contours  of  the  global  error  as  a  function  of  op  and  To-  The  values  of  |j-  between  1 
to  2  do  not  significantly  affect  the  overall  accuracy.  The  optimal  value  of  To  is  around  10  ~  12. 

(b)  Global  etTor  variation  with  respect  to  To-  PML  thickness  N  and  S  velocities  cs  (ao  is  fixed 
to  7 r/c).  The  mimimal  global  error  is  around  11  for  different  PML  thicknesses  of  10,  12  and  20 
(c.s  =  1600  m/s),  indicating  the  optimal  To  is  not  sensitive  to  the  thickness  of  the  PML.  For  N  =  20 
and  cs  =  2200  m/s,  the  optimal  To  is  around  15.  (c)  Normalized  points  per  dominant  wavelength 
variation  inside  the  PML  layer  for  S  wave  for  different  PML  thickness  N  and  S  wave  speed  cs.  (d) 
The  zoom  view  of  Figure  8c  to  highlight  details.  The  normalized  points  per  dominant  wavelength 
is  around  0.5  at  the  exterior  boundary  in  difference  cases. 

9  Thin  slab  model  with  the  free  surface.  A  10-cell-thick  PML  (dashed  lines)  is  applied 
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on  all  the  faces  except  the  free  surface.  The  star  symbol  shows  the  location  of  the  vertical  force 
(200,  3000,  795)  m,  ten  cells  away  from  the  PML-interior  interface  at  the  negative  x  direction  and 
half  cell  below  the  free  surface.  Other  parameters  are  the  same  as  in  Figure  3. 

10  Snapshots  of  the  vz  component  of  the  free-surface  slab  model  at  3.75  s.  The  color  is  scaled 
to  0.1%  of  the  maximum  absolute  amplitude.  White  lines  show  the  PML-interior  interfaces,  (a)  Ref¬ 
erence  solution,  (b)  Standard  PML.  (c)  /5o  =  1  and  ao  =  irfc.  (d)  Standard  PML  without  using  the 
modified  PML  equations  31  and  32.  Strong  spurious  reflected  surface  waves  can  be  clearly  seen. 

(e)  Standard  PML  with  the  free  surface  coincident  with  the  location  of  the  crxz  component. 

11  The  two-layered  model  with  the  free  surface  topography,  which  is  defined  by  z  =  1000  exp(— r2/10002), 
where  r  =  \J x'2  +  y2.  When  r  >  2000m,  zq  is  damped  to  zero  by  a  Gaussian  function.  The  model 

is  19.95km  in  length,  6.15km  in  width  and  10.7km  in  depth  (including  PML  layers).  The  interface 
is  at  z  =  —8km,  and  the  medium  parameters  are  cp  =  3000  m/s,  cs  =  2000  m/s,  p  =  1200  kg/m3 
for  the  upper  layer  and  cp  =  5000  m/s,  cs  =  2500  m/s,  p  =  1800  kg/m3  for  the  lower  layer.  A 
12-cell-thick  PML  is  applied  on  all  the  faces  except  the  free  surface  (white  lines).  The  star  symbol 
shows  the  location  of  the  vertical  force  (—8000, 0,  —500)  m,  and  the  triangle  symbols  indicates  the 
receiver’s  location  (9000,  0,  0)  m,  seven  cells  distance  to  the  PML-interior  interface  at  the  positive 
x  face. 

12  Snapshots  of  vz  at  7  s  on  the  xz  cross-section  throuhg  the  middle  point  of  y  axis.  The  color 
is  scaled  to  less  than  1%  of  the  maximum  absolute  amplitude.  White  lines  are  the  PML-interior  in¬ 
terfaces.  (a)  Reference  solution,  (b)  (3q  =  1,  «o  =  vr/c. 

13  (a)  Comparison  of  seismograms  (vz)  for  the  reference  (solid  line),  So  =  1  (dashed  line) 
and  Pq  =  8  (dotted  line)  at  the  receiver  shown  in  Figure  11.  (b)  Detailed  comparison  within  the  time 
window  13  ~  21  s  and  the  amplitude  is  zoomed  in  to  0.2%  of  that  in  Figure  13a.  Difference  can 
be  seen  between  the  reference  and  /5q  =  1.  The  solution  with  So  =  8  overlaps  with  the  reference. 
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(c)  Difference  between  solutions  from  the  CFS  ADE-PML  and  the  reference.  Solid  line  shows  the 
difference  for  (3q  =  1,  and  dotted  line  for  the  optimal  /3q  =  8. 
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Figure  1:  A  3D  computational  domain  surrounded  by  PML  layers  on  each  face.  Dashed  lines  mark 
the  PML-interior  interfaces.  Along  the  edge  of  the  computational  domain,  two  directional  PMLs 
exist.  At  the  corner,  three  directional  PMLs  co-exist. 
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(a)  (b) 


Figure  2:  (a)  Maximum  amplitude  distribution  for  different  (3  values  (/?o  =  solid  line  and  /3o  =  9, 
dotted  line)  for  a  25 -cell-thick  PML  of  a  normal  incident  plane  wave,  a  =  0  and  an  optimal  d 
are  used,  (b)  Velocity  seismogram  at  the  grid  point  15-cell  into  the  PML  for  /?o=  1  (solid  line)  and 
/3q  =  9  (dotted  line).  The  wave  using  /3q  =  9  arrives  later  than  that  using  Bo  =  1  due  to  the  reduced 
wave  speed  by  /3X- 
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Figure  3:  A  thin  homogeneous  3D  slab  surrounded  by  10-cell-thick  PMLs  on  all  six  faces  (three 
cross-sections  are  shown).  The  PML-interior  interfaces  are  indicated  by  dashed  lines.  The  medium 
has  a  compressional  wave  speed  cp  =  3200  m/s,  shear  wave  speed  cs  =  1600  m/s,  and  density 
p  =  2200  kg/m3.  The  grid  spacing  is  10  m,  the  total  grid  number  is  600  x  600  x  81,  and  the  size 
of  the  model  is  5990  x  5990  x  800  m.  The  star  symbol  shows  the  location  of  the  vertical  force 
(200,  3000,  600)  m,  ten  cells  away  both  from  the  PML-interior  interface  at  the  negative  x  face  and 
from  that  at  the  positive  z  face.  The  triangle  symbol  indicates  the  receiver  (5790,  3000,  600)  m,  ten 
cells  away  from  both  the  PML-interior  interface  at  the  positive  x  face  and  that  at  the  positive  z  face. 
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(e) 


(0 


Figure  4:  Snapshots  of  the  v~  component  of  the  thin  slab  model  (Figure  3)  at  3.75  s.  The  color  is 
scaled  to  0.1%  of  the  maximum  absolute  amplitude.  White  lines  are  the  PML-interior  interfaces, 
(a)  Reference  solution,  (b)  Standard  PML.  (c)  oq  =  7 r/c  and  ;3q  =  1.  (d)  o:o  =  7 r/c  and  /?o  =  10. 
(e)  ao  =  7 r/c,  f3o  =  7  and  a  three  times  larger  do-  (f)  ao  =  vr/c  and  /3o  =  11  for  a  20-cell-thick 
PML. 
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Figure  5:  The  vz  seismograms  at  the  receiver  shown  in  Figure  3.  The  right  column  shows  the 
same  seismograms  except  that  the  y  axis  is  scaled  to  1%  of  that  of  the  left  column  to  highlight 
the  difference.  The  first  arrival  in  the  right  column  is  P  wave,  which  is  nearly  invisible  in  the  left 
column,  (a)-(b)  Standard  PML.  (c)-(d)  /?o  =  1  and  ao  =  vr  fc.  (e)-(f)  /3q  =  10  and  ao  =  i r/c.  (g)-(h) 
/?o  =  7,  oq  =  7r  fc  and  a  three  times  larger  do-  41 
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x-axis  (m)  x-axis  (m) 


(a) 
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Figure  6:  Snapshots  of  vz  on  the  cross-section  in  the  middle  of  y  axis  at  3.75  s.  The  color  is 
scaled  to  the  same  maximum  absolute  amplitude.  Green  lines  are  the  PML-interior  interfaces,  (a) 
a0  =  vr/c  and  /3o  =  1.  The  amplitude  of  S  wave  in  the  upper  PML  layer  is  even  larger  than  that 
in  the  physical  region,  which  indicates  a  does  not  efficiently  absorb  this  near-grazing  propagating 
body  wave.  The  distribution  of  the  peak  amplitude  is  wider  than  Figure  6b,  which  indicates  the 
reflected  wave  propagates  into  the  physical  domain,  (b)  qq  =  7 r/c  and  (3o  =  10.  The  wave  front  is 
bended  towards  the  normal  direction  inside  both  the  upper  PML  layer  and  the  bottom  PML  layer. 
The  amplitude  of  S  wave  exhibits  damping  in  the  PML  layers. 
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(C) 


(d) 


Figure  7:  Local  error  distribution.  The  color  is  scaled  to  0.5.  White  lines  show  the  PML-interior 
interfaces,  (a)  Standard  PML.  (b)  /3q  =  1  and  «o  =  tt/c-  (c)  (3q  =  10  and  ao  =  7 r/c.  (d)  3q  =  7, 
q:(i  =  7 r/c  and  a  three  times  larger  do- 
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(a)  (b) 


(c)  (d) 


Figure  8:  (a)  Contours  of  the  global  error  as  a  function  of  oo  and  /3q.  The  values  of  “j.1  between 
1  to  2  do  not  significantly  affect  the  overall  accuracy.  The  optimal  value  of  /3q  is  around  10  ~  12. 
(b)  Global  error  variation  with  respect  to  (3q ,  PML  thickness  N  and  S  velocities  cs  (on  is  fixed 
to  7T  fc).  The  mimimal  global  error  is  around  11  for  different  PML  thicknesses  of  10,  12  and  20 
(cs  =  1600  m/s),  indicating  the  optimal  (3q  is  not  sensitive  to  the  thickness  of  the  PML.  For  N  =  20 
and  cs  =  2200  m/s,  the  optimal  3q  is  around  15.  (c)  Normalized  points  per  dominant  wavelength 
variation  inside  the  PML  layer  for  S  wave  for  different  PML  thickness  N  and  S  wave  speed  cs.  (d) 
The  zoom  view  of  Figure  8c  to  highlight  details.  The  normalized  points  per  dominant  wavelength 
is  around  0.5  at  the  exterior  boundary  in  difference  cases. 
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Figure  9:  Thin  slab  model  with  the  free  surface.  A  10-cell-thick  PML  (dashed  lines)  is  applied 
on  all  the  faces  except  the  free  surface.  The  star  symbol  shows  the  location  of  the  vertical  force 
(200,  3000,  795)  m,  ten  cells  away  from  the  PML-interior  interface  at  the  negative  x  direction  and 
half  cell  below  the  free  surface.  Other  parameters  are  the  same  as  in  Figure  3. 

Zhang  &  Shen  - 


45 


A-46 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 

21 

22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 


GEOPHYSICS 


Page  46  of  49 


(e) 


Figure  10:  Snapshots  of  the  vz  component  of  the  free-surface  slab  model  at  3.75  s.  The  color  is 
scaled  to  0.1%  of  the  maximum  absolute  amplitude.  White  lines  show  the  PML-interior  interfaces. 


(a)  Reference  solution,  (b)  Standard  PML.  (c)  [3q  =  1  and  cto  =  7r/c.  (d)  Standard  PML  without 
using  the  modified  PML  equations  3 1  and  32.  Strong  spurious  reflected  surface  waves  can  be  clearly 
seen,  (e)  Standard  PML  with  the  free  surface  coincident  with  the  location  of  the  axz  component. 
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Figure  1 1 :  The  two-layered  model  with  the  free  surface  topography,  which  is  defined  by  z  = 
1000  exp(— r2/10002),  where  r  =  \J  x2  +  y2.  When  r  >  2000m,  zq  is  damped  to  zero  by  a  Gaus¬ 
sian  function.  The  model  is  19.95km  in  length,  6.15km  in  width  and  10.7km  in  depth  (including 
PML  layers).  The  interface  is  at  z  =  —8km,  and  the  medium  parameters  are  cp  =  3000  m/s, 
cs  =  2000  m/s,  p  =  1200  kg/m3  for  the  upper  layer  and  cp  =  5000  m/s,  cs  =  2500  m/s, 
p  =  1800  kg/m3  for  the  lower  layer.  A  12-cell-thick  PML  is  applied  on  all  the  faces  except  the  free 
surface  (white  lines).  The  star  symbol  shows  the  location  of  the  vertical  force  (—8000, 0,  —500)  m, 
and  the  triangle  symbols  indicates  the  receiver’s  location  (9000, 0, 0)  m,  seven  cells  distance  to  the 
PML-interior  interface  at  the  positive  x  face. 
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Figure  12:  Snapshots  of  vz  at  7  s  on  the  xz  cross-section  throuhg  the  middle  point  of  y  axis.  The 
color  is  scaled  to  less  than  1%  of  the  maximum  absolute  amplitude.  White  lines  are  the  PML- 
interior  interfaces,  (a)  Reference  solution,  (b)  /3q  =  1,  op  =  tt/c- 
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Figure  13:  (a)  Comparison  of  seismograms  (vz)  for  the  reference  (solid  line),  0o  =  1  (dashed  line) 
and  I30  =  8  (dotted  line)  at  the  receiver  shown  in  Figure  11.  (b)  Detailed  comparison  within  the  time 
window  13  ~  21  s  and  the  amplitude  is  zoomed  in  to  0.2%  of  that  in  Figure  13a.  Difference  can 
be  seen  between  the  reference  and  0q  =  1.  The  solution  with  So  =  8  overlaps  with  the  reference, 
(c)  Difference  between  solutions  from  the  CFS  ADE-PML  and  the  reference.  Solid  line  shows  the 
difference  for  0q  =  1,  and  dotted  line  for  the  opt^fial  So  =  8. 
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SUMMARY 

Advance  in  computation  has  made  routine  use  of  Green  functions  based  on  three-dimensional 
(3D)  reference  models  an  increasingly  attractive  option  in  earthquake  source  inversion,  as  3D 
reference  models  help  to  isolate  source  phenomena  from  propagation  effects.  This  paper  extends 
the  source  inversion  method  based  on  the  strain  Green  tensor  (SGT)  database  from  local  to 
regional  scales  and  presents  the  location  and  moment-tensor  determinations  of  earthquakes 
between  2001  and  2004  in  the  southeast  Tibetan  plateau.  A  non- staggered  grid  finite-difference 
(FD)  method  is  used  to  simulate  wave  propagation  in  a  3D  spherical  earth  model  and  construct 
the  SGT  database.  The  FDSGT  inversion  takes  advantage  of  the  source-receiver  reciprocity  and 
is  well  suited  for  (near)  real  time  earthquake  source  studies.  The  SGT  database  approach  makes 
it  computationally  efficient  to  determine  the  best  solution  in  a  grid  search  surrounding  the 
reference  location,  by  minimizing  the  waveform  misfit  of  time-windowed  first  arrivals, 
teleseismic  S  waves  and  surface  waves  recorded  at  local  to  teleseismic  distances.  The 
waveforms  are  bandpass-filtered  in  four  frequency  bands  to  accommodate  frequency-dependent 
waveform  quality  and  time  shifted  to  accommodate  errors  in  the  SGTs  due  to  an  imperfect 
reference  model.  No  prior  source  mechanism  solution  is  needed  in  the  inversion.  The  new 
solutions  are  in  general  agreement  with  the  global  centroid  moment  tensor  solutions,  though  the 
two  source  locations  may  differ  up  to  10s  km.  The  moment  tensor  solutions  are  not  sensitive  to 
uncertainties  in  lateral  source  location,  partly  due  to  the  time  shift  of  waveforms.  Inspection  of 
waveform  fit  in  detail  suggests  that  regional  stations,  the  3D  reference  model  and  short-period 
body  waves  tend  to  provide  tighter  constraints  on  the  locations  of  the  earthquakes.  The  lowest 
magnitude  of  the  earthquakes  in  this  study  is  Mw  4.35.  As  the  crustal  velocity  model  improves, 
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one  may  extend  useful  waveforms  to  short-period  (<15  s)  surface  waves  and  be  able  to  apply  the 
FDSGT  method  to  smaller  earthquakes  at  regional  scales. 

Key  words:  moment  tensor  inversion,  earthquake  source  mechanism,  Tibetan  plateau,  three- 
dimensional  reference  models,  strain  Green  tensor  database 

1  INTRODUCTION 

Accurate  source  mechanisms  are  essential  in  the  studies  of  earthquakes  as  well  as  waveform- 
based  investigations  of  the  earth  structure.  Earthquake  focal  mechanisms  can  be  determined 
from  first-motion  polarities  or  a  combination  of  first-motion  polarities  and  ratios  of  maximum  P 
amplitude  to  maximum  S  amplitude  (e.g.,  Kisslinger,  1980;  Shen  et  al.,  1997;  Hardebeck  & 
Shearer,  2003).  With  the  development  of  digital  seismograph,  seismologists  have  also  pursued 
various  waveform  approaches  in  source  inversion  by  minimizing  the  misfit  between  synthetic 
and  observed  body  and  surface  waveforms  (e.g.,  Dziewonski  et  al.,  1981;  Nabelek,  1984;  Dreger 
&  Helmberger,  1991;  Kikuchi  &  Kanamori,  1991;  Zhao  &  Helmberger,  1994;  Zhu  & 
Helmberger,  1996).  While  source  mechanisms  of  moderate-to-large  earthquakes  have  been 
routinely  determined  by  modeling  teleseismic  waveforms  (e.g.,  Dziewonski  et  al.,  1981), 
characterization  of  small-to-moderate  earthquakes  has  been  uneven  geographically  and  usually 
incomplete.  Yet  small-to-moderate  earthquakes  provide  valuable  information  about  regional 
tectonics  due  to  their  relative  abundance  and  broader  distribution.  Part  of  the  difficulties  in 
applying  a  generic  moment  tensor  inversion  scheme  to  different  areas  is  the  strong  lateral 
variation  in  crustal  structure,  which  generates  complex  short-period  waves  propagating  at  local 
and  regional  distances. 
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Most  earthquake  source  studies  to  date  utilize  one-dimensional  (ID)  earth  model  in  the 
calculation  of  synthetic  waveforms  or  Green  functions.  To  accommodate  three-dimensional 
(3D)  structural  heterogeneities,  segments  of  body  and  surface  waves  are  time-shifted  before 
waveform  fitting  (e.g.,  Zhao  &  Helmberger,  1994;  Zhu  &  Helmberger,  1996).  In  places  with 
highly  heterogeneous  earth  structure,  this  “cut-and-paste”  of  waveform  segments  may  not  fully 
account  for  the  complexity  of  waveform  phenomena  caused  by  3D  structures  (e.g.,  focusing  and 
defocusing  of  waves).  One  of  the  major  sources  of  errors  in  moment  tensor  inversion  is  phase 
skipping  between  the  observed  and  synthetic  waveforms.  When  traveltime  delays  due  to  3D 
velocity  heterogeneities  are  comparable  to  the  wave  period,  phase  skipping  becomes  a 
challenging  problem. 

Advance  in  computation  has  made  a  routine  use  of  3D  Green  function  in  earthquake  source 
inversion  an  increasingly  attractive  option.  Synthetic  waveforms  from  good  3D  reference 
models  account  for  complex  waveform  phenomena  and  reduce  phase  differences.  Liu  et  al. 
(2004)  determined  the  source  mechanisms  of  small  to  moderate  earthquakes  in  southern 
California  using  synthetic  waveforms  calculated  with  the  spectral-element  method.  Because  the 
derivatives  of  the  source  parameters  are  determined  numerically  by  differentiating  synthetics 
with  respect  to  the  source  parameters,  up  to  10  forward  simulations  are  needed  for  each 
earthquake  (6  moment  tensor  components,  latitude,  longitude,  depth,  and  the  reference  location). 
When  the  reference  location  is  far  from  the  true  location,  the  derivatives  of  the  location 
parameters  may  not  adequately  reflect  the  non-linear  variation  of  waveforms  as  a  function  of 
location.  Taking  advantage  of  source-receiver  reciprocity,  Zhao  et  al.  (2006)  introduced  the  use 
of  strain  Green  tensors  (SGTs)  in  source  inversion  based  on  3D  reference  models.  In  this 
approach,  the  SGTs  are  calculated  with  a  finite-difference  code  and  saved  in  a  database.  The 


B-5 


4 


Page  5  of  42 


Geophysical  Journal  International 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 
11 
12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 

43 

44 

45 

46 

47 

48 

49 

50 

51 

52 

53 

54 

55 

56 

57 

58 

59 

60 


number  of  forward  waveform  simulations  is  three  times  the  number  of  stations  used  in  source 
inversion.  When  the  sources  outnumber  the  stations,  the  Green’s  tensor  or  SGT-database 
approach  is  more  computationally  efficient  (e.g.,  Graves  &  Wald,  2001;  Eisner  &  Clayton, 
2001).  More  importantly  since  the  SGTs  are  pre-calculated,  the  SGT-database  approach  is  better 
suited  for  (near)  real-time  source  inversion.  Zhao  et  al.  (2006)  adopted  a  perturbation  approach 
in  source  inversion  and  chose  to  minimize  frequency-dependent  arrival  time  and  amplitude 
anomalies. 

In  this  paper,  we  extend  the  SGT-database  approach  from  local  scales  with  waveforms 
calculated  in  a  Cartesian  coordinate  (Zhao  et  al.,  2006)  to  regional  and  teleseismic  scales  with 
waveforms  calculated  in  a  spherical  coordinate  (Zhang  et  al.,  in  preparation),  so  we  are  able  to 
use  any  combination  of  local  and  teleseismic  waveforms  in  source  inversion.  Instead  of 
frequency-dependent  arrival  time  and  amplitude  anomalies,  we  fit  multi-frequency  filtered  and 
time-shifted  waveforms  directly.  For  a  given  source  location,  the  waveforms  are  linear  functions 
of  the  moment  tensor  elements,  which  can  be  solved  as  a  linear  inverse  problem.  To  account  for 
perturbations  in  earthquake  source  location,  we  carry  out  a  grid  search  in  the  vicinity  of  the 
reference  location.  This  is  necessary  for  earthquakes  with  location  errors  of  5  km  and  more,  as 
in  the  case  of  most  global  earthquake  catalogs.  When  the  true  source  location  is  far  from  the 
reference  location,  the  derivative  of  waveform  misfit  about  the  reference  location  may  not  be  a 
good  approximation  to  obtain  the  true  location.  This  in  turn  affects  the  waveforms  and  the 
moment  tensor  solution.  The  grid  search  is  computationally  efficient  as  it  involves  extracting 
pre-calculated  SGTs  in  a  small  volume  containing  the  reference  source  location  from  the 
database.  We  apply  the  new  method  to  southeast  Tibet  in  conjunction  with  a  full-wave 
tomographic  imaging  effort,  as  good  earthquake  source  mechanisms  are  essential  for  waveform- 
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based  tomographic  inversions.  An  updated  upper  mantle  model  (CUB  2.0,  Shapiro  & 
Ritzwoller,  2002)  is  used  as  the  3D  reference  model. 

2  SEISMOTECTONICS  OF  THE  TIBETAN  PLATEAU  AND  SEISMIC  DATA 

Due  to  the  collision  between  the  Indian  and  Eurasian  plates,  the  Tibetan  plateau  has  abundant 
and  widely  distributed  seismic  activities.  Source  mechanism  solutions  show  diverse  styles  of 
thrust,  normal  and  strike-slip  faulting  in  the  plateau  and  surrounding  areas  (e.g.,  Chen  et  al., 
1981;  Zhou  et  al.,  1983;  Ni  &  Barazangi,  1984;  Pettersen  &  Doornbos,  1987;  Wu  et  al.,  2004). 
Using  seismic  waves  recorded  by  the  WWSSN  seismographs,  Monlar  &  Lyon-Caen  (1989)  and 
Zheng  &  Suzuki  (1993)  studied  the  moment  tensors  of  earthquakes  occurred  in  the  region  in  the 
60s  to  80s.  Randall  et  al.  (1995)  used  regional  seismograms  recorded  during  the  1991-1992 
PASSCAL  experiment  on  the  Tibetan  plateau  and  estimated  the  source  mechanisms  for 
earthquakes  with  Mw  ranging  from  3.8  to  5.0.  Similarly  regional  earthquake  records  from  other 
PASSCAL  experiments  in  central  Tibetan  plateau  (Langin  et  al.,  2003)  and  southern  Tibet  (de  la 
Torre  et  al.,  2007)  were  used  to  determine  earthquake  source  mechanisms.  In  this  study,  we 
focus  on  earthquakes  in  the  southeastern  Tibetan  plateau  recorded  by  the  HIMNT  (de  la  Torre  & 
Sheehan,  2005),  BUTHAN  (Velasco  et  al.,  2007),  Namche  Barwa  (Sol  et  al.,  2007),  and  MIT- 
China  (Xu  et  al.,  2007)  experiments  between  2001  and  2004.  We  augment  seismic  data  recorded 
by  portable  instruments  with  records  from  selected  permanent  seismic  stations  at  regional  and 
teleseismic  distances  (Fig.  1). 

3  MOMENT  TENSOR  INVERSION 
3.1  The  inverse  problem 
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Following  Kikuchi  &  Kanamori  (1991),  we  represent  a  general  moment  tensor  M  as  a  linear 
combination  of  6  elementary  (basis)  moment  tensors, 

6 

m=v«a  (i) 

mm  v  ' 

m=  1 

where  am  is  the  coefficient  for  the  elementary  moment  tensors  and, 
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We  adopt  this  representation  because  of  an  important  advantage  of  using  the  elementary 
moment  tensors  instead  of  the  6  independent  elements  of  a  general  moment  tensor.  Namely, 
specific  solutions  (e.g.,  double  couple  with  a  vertical  nodal  plane  or  pure  deviatoric  moment 
tensor)  can  be  obtained  with  subgroups  of  the  elementary  moment  tensors,  thus  requiring  fewer 
than  6  coefficients  or  unknowns  in  the  inversion.  Let  b]nm{t,p)  denote  the  displacement 
component  n  at  station  j  from  an  elementary  moment  tensor  Mm  at  source  location  p.  b]nm(t,p) 
can  be  constructed  from  the  SGTs  for  the  station  (Zhao  et  al.,  2006), 

bUt,p)=jj.Gu(p,r,j)M*]tS(t,p) 

i,k=i  (2) 

i,k= 1 
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where  E"k(p,t;j)  is  the  strain  tensor  component  ik  at  source  location  p  due  to  a  unit  impulsive 
force  acting  at  receiver  j  in  the  direction  n,  M*  the  component  ik  of  the  elementary  moment  Mm, 
and  S(t,  p)  the  source  time  function. 

The  displacement  s„  due  to  a  general  moment  tensor  can  be  expressed  as  a  linear  combination 
of  the  basis  displacements, 

sJn(t,p)  =  iccmb/im(t,p)  (3) 

m=  1 

where  Nb  is  the  number  of  the  elementary  moment  tensors  used  (up  to  6).  For  brevity,  we 
hereinafter  drop  the  superscript  denoting  receiver  and  the  subscript  denoting  the  direction  of  the 
displacement  and  define  a  least-square  misfit  function  as, 

E(  m,  p)  =  wt  Jj  K(0  -  iambim(E  p)fdt  (4) 

i=l  m=  1 

where  J,(t)  is  the  z-th  seismogram,  N  the  number  of  time  segments,  ti  to  h  the  time  window  of 
the  waveform,  m  =  \av...,aN  ]rthe  unknown  coefficients  to  be  solved,  and  vv/  a  weighting  factor 
to  account  for  waveform  quality,  station  azimuthal  distribution  and  amplitude.  The  misfit 

ZN 

wr  The  normal  equation  for  the  least-square  problem  is 

(Kikuchi  &  Kanamori,  1991), 

Nb 

EcA  =  D„  for  k=l,...,Nb  (5) 

m=  1 

where 


Gkm  =  j  WE  P)bim(t,  p)\dt 


i=\ 


(6) 
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D,  = 


N  f 

Yjw,  J  [bik(t,p)di(t)]dt 


(7) 


The  solution  m  =  [av...,aNb  ]T  is  obtained  by  singular  value  decomposition  and  the  moment 
tensor  solution  is  calculated  from  equation  1. 


3.2  The  3D  reference  model  and  construction  of  strain  Green’s  tensor  database 

One  way  to  select  a  3D  reference  model  is  to  compare  the  goodness  of  fit  between  the 
observed  and  synthetic  waveforms.  We  compare  CRUST  2.0  (Bassin  et  al.,  2000)  and  CUB  2.0 
(Shapiro  &  Ritzwoller,  2002)  in  this  study.  We  extend  CRUST  2.0  beneath  the  Moho  with  the 
AK135  model  (Kennett  et  al.,  1995)  and  scale  P-wave  speed  to  S-wave  speed  in  CUB  2.0  using  a 
Vp/Vs  of  1.74  in  the  crust  and  that  of  AK135  in  the  mantle.  The  density  is  converted  from  the 
density-velocity  relationship  of  Christensen  &  Mooney  (1995).  When  necessary,  we  also  extend 
CUB  2.0  to  below  400  km  depth  with  AK135.  At  the  reference  model  selection  stage,  we  use 
the  global  centroid  moment  tensor  solutions  (GCMT,  formerly  known  as  Harvard  CMT)  of  select 
earthquakes  to  calculate  synthetic  waveforms.  The  band-pass  filtered  synthetic  waveforms  of  the 
entire  earthquake  wave  train  are  cross-correlated  with  the  observed.  Fig.  2  shows  that  on 
average  CRUST  2.0  predicts  a  much  earlier  wave  train  than  the  observed  and  more  importantly  a 
poor  maximum  cross-correlation  coefficient  than  CUB  2.0  at  all  frequency  bands.  As  expected, 
the  correlation  coefficients  are  close  to  1  and  relatively  less  scattered  at  long  periods  (50-120  s). 
At  periods  of  15  s  and  shorter,  the  correlation  coefficients  are  highly  variable.  Since  the 
correlation  is  dominated  by  surface  waves,  which  have  relatively  large  amplitudes  compared  to 
those  of  body  waves,  the  highly  variable  correlation  coefficients  at  short  periods  indicates  that 
these  two  models  do  not  accurately  represent  the  real  shallow  crust. 
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We  select  CUB  2.0  as  our  initial  reference  model  in  source  mechanism  and  full-wave 
tomographic  inversions  for  its  overall  better  performance  in  the  waveform  cross-correlation  test 
(Shen  et  al.,  2008).  The  source  mechanism  inversion  presented  in  this  paper  uses  a  new  model 
derived  from  CUB  2.0  after  the  first  iteration  of  full-wave  inversion. 

To  construct  the  station  SGTs  to  be  used  in  source  inversion,  we  use  a  non-staggered-grid 
finite-difference  method  in  the  spherical  coordinate  (Zhang  et  al.,  in  preparation).  Three 
simulations,  one  for  each  orthogonal  unit  impulsive  point  force  acting  at  the  receiver,  are  needed 
for  each  station  used  in  the  source  inversion.  We  select  49  stations  at  regional  and  teleseismic 
distances  from  the  temporary  and  permanent  seismic  networks  (Fig.  1).  Simulations  are  carried 
out  on  three  geographic  scales:  (1)  within  the  southeast  Tibetan  plateau,  (2)  at  far-regional 
distances,  and  (3)  teleseismic  distances.  The  finite-difference  grid  in  the  regional  and  far- 
regional  simulations  has  a  grid  spacing  of  0.04°  (latitude),  0.045°  (longitude),  and  1-6  km  (with 
the  small  depth  interval  near  the  free  surface).  This  yields  accurate  P  wave  synthetics  at  periods 
great  than  ~7  s  and  accurate  shear  and  surface  wave  synthetics  at  periods  greater  than  -10  s.  The 
moment  rate  is  a  Gaussian  time  function  with  an  effective  cutoff  frequency  below  ~0.5  Hz  and 
centered  at  3  s.  Because  the  synthetics  are  low-pass  filtered  to  -0.13  Hz,  the  periods  of  the 
waveforms  used  in  the  inversion  are  much  longer  than  the  typical  source  duration  of  the  earthquakes 
in  this  study  (M  4-6).  For  larger  earthquakes,  it  may  be  necessary  to  use  an  empirically  determined 
source  time  magnitude  relationship  or  determine  the  source-time  function  in  the  source  inversion. 
For  each  simulation,  we  use  16  cpu  cores  (AMD  Opteron  2218,  2.6  Ghz,  2  GB  memory  per  core) 
and  carry  out  up  to  6  simultaneous  simulations  on  a  104-core  cluster.  The  finite-difference  grid 
in  the  teleseismic  simulations  has  three  times  the  lateral  spacing  and  twice  the  vertical  spacing  of 
the  regional  simulations.  Thus  for  teleseismic  arrivals,  we  use  only  synthetics  at  periods  greater 
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than  30  s.  The  SGTs  are  saved  on  the  same  spatial  grids  (0.12°  in  latitude,  0.135°  in  longitude, 
and  2  km  in  depth)  within  the  southeast  Tibetan  plateau,  because  all  earthquakes  studied  in  this 
paper  are  within  the  southeast  Tibetan  plateau.  The  same  SGTs  are  also  used  to  calculate  the 
finite-frequency  sensitivity  kernels  (Zhao  et  al.,  2005;  Zhang  et  al.,  2007)  in  a  parallel  effort  to 
improve  the  velocity  model  by  full-wave  tomography  inversion.  For  a  study  just  on  earthquake 
moment  tensors,  one  can  save  SGTs  only  in  the  vicinity  of  the  target  earthquakes  and  in  finer 
grids.  The  SGTs  extracted  from  the  database  for  the  inversions  of  the  23  earthquakes  in  this 
study  take  ~23  GB  disk  space. 

3.3  Data  Processing 

From  a  list  of  earthquakes  in  the  southern  Tibetan  plateau  between  2001  and  2004  in  the 
EHB  earthquake  catalog  (E.R.  Engdahl,  personal  communication),  we  collect  earthquake  records 
from  the  portable  and  permanent  broadband  seismic  stations  in  southeastern  Tibet  and  selected 
stations  at  far-regional  and  teleseismic  distances  (Fig.  1).  We  remove  the  instrument  response 
and  bandpass  filter  the  records  in  four  frequency  bands  (0.00833-0.0167,  0.0167-0.0333,  0.0333- 
0.0667,  0.0667-0.133  Hz).  We  select  the  records  by  their  signal-to-noise  ratios  (SNRs)  using  an 
automated  process.  The  signal  level  is  the  maximum  amplitude  of  the  trace,  while  the  noise  level 
is  defined  as  the  amplitude  of  the  trace  in  the  time  window  between  100  s  prior  to  the  earthquake 
origin  time  and  the  first  arrival.  Those  with  SNRs  greater  than  a  cut-off  value  (3)  are  selected  for 
further  analysis.  Bandpass-filtering  allows  records  useable  in  one  frequency  band  but  noisy  in 
other  frequency  bands  to  pass  the  SNR  test.  We  further  select  the  earthquakes  by  their 
geographic  distributions  and  limit  our  source  inversions  to  earthquakes  that  are  well  recorded  by 
at  least  5  stations. 

3.3.1  Selection  of  waveforms 
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Following  Zhao  &  Helmberger  (1994),  we  select  two  time  windows  on  each  seismic  record. 
At  regional  and  far-regional  distances,  we  use  the  surface  wave  and  the  first  arrival  (Pg,  Pnl), 
though  we  avoid  the  first  arrival  at  the  upper  mantle  triplication  distance  (1400-3000  km).  This 
restriction  can  be  relaxed  if  seismic  records  at  other  distances  are  insufficient  and  the  benefit  of 
the  addition  of  body  waves  out-weights  the  potential  downside  of  an  inaccurately  modeled 
complex  triplication  waveform.  At  teleseismic  distances  (>  3000  km),  we  use  the  first  arrival  (P) 
and  S  wave.  Generally  the  first  arrivals  and  teleseismic  body  waves  are  less  affected  by  shallow 
heterogeneities  and  provide  relatively  consistent  waveforms  between  the  synthetics  and 
observed,  though  they  tend  to  have  lower  SNRs  due  to  their  relatively  low  amplitude.  Long- 
period  (60-120  s)  surface  waves  also  provide  consistent  waveforms.  But  for  smaller  earthquakes 
(M<5),  intermediate-period  records  often  have  better  SNRs. 

Because  synthetic  surface  waves  from  the  current  reference  model  do  not  fit  the  observed 
well  at  short  periods  (15  s  and  less,  Fig.  2),  we  limit  surface  waves  to  30  s  and  longer  period  in 
the  source  inversion.  As  the  velocity  model  improves,  one  can  increase  the  bandwidth  of  useable 
surface  waves  to  shorter  periods.  We  also  exclude  the  lowest  frequency  body  waves  (since  they 
often  overlap  with  surface  waves  at  long  period)  and  only  use  body  wave  signals  at  periods 
below  60  s.  For  the  first  arrival,  the  length  of  the  time  window  varies  with  the  frequency  band. 
From  the  onset  of  the  first  arrival  to  the  end  of  the  time  window  is  one  full  period  of  the  long- 
period  end  of  the  frequency  band.  This  is  designed  to  minimize  contributions  from  later  arrivals 
(e.g.,  PmP),  which  are  likely  less  well  captured  by  the  reference  model. 

3.3.2  Waveform  weighting  and  time  shift 

The  3D  reference  model  helps  to  reduce  the  phase  shifts  between  the  observed  and 
synthetics,  but  significant  phase  mismatches  remain  due  to  errors  and  unaccounted 
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heterogeneities  in  the  reference  model.  As  in  Zhao  &  Helmberge  (1994),  we  allow  the  observed 
and  synthetic  waveforms  to  shift  in  time  and  calculate  the  misfit  between  the  observed  and  time 
shifted  synthetics.  The  amount  of  time  shift  is  determined  by  waveform  cross-correlation  and 
limited  within  a  fraction  of  the  arrival  time  of  the  phase.  This  is  equivalent  to  a  prescribed  value 
that  varies  with  distance  (Liu  et  al.,  2004).  By  trial  and  error,  we  set  the  allowable  time  shift  of  a 
body  wave  to  within  ±3%  of  its  travel  time.  For  surface  waves,  we  allow  the  phase  to  shift  by  a 
larger  fraction  (±5%).  Too  large  a  fraction  may  lead  to  phase  skip  and  unreasonably  large  phase 
shifts,  while  too  small  a  fraction  does  not  allow  otherwise  very  similar  waveforms  to  match  in 
time.  As  the  velocity  model  improves,  one  may  use  smaller  fractional  values,  which  are 
expected  to  generate  a  stronger  gradient  in  the  misfit  function. 

Since  data  quality  and  source-receiver  distribution  vary  greatly  from  earthquake  to 
earthquake,  from  station  to  station,  and  from  frequency  band  to  frequency  band,  it  is  often 
necessary  to  treat  the  data  by  a  weighting  scheme  that  takes  into  account  our  knowledge  of 
individual  data  quality.  Different  researchers  may  adopt  different  weighting  schemes,  making 
this  perhaps  the  most  subjective  aspect  of  source  inversion.  Though  for  earthquakes  with  good 
azimuthal  coverage  and  high  SNR  records,  weighting  is  likely  non-critical. 

The  assigned  weight  in  equation  4,  vv,  =  w'"r  w',zw'i‘",p ,  is  composed  of  weighting  for  signal 
noise  ratio  w-nr ,  station  distribution  vv"" ,  and  amplitude  w“mp .  To  account  for  various  qualities  of 
seismic  arrivals,  we  assign  w-nr= 1  for  time  windows  with  SNR  between  3-6  and  vv”'r=2  for 
windows  with  SNR  >  6.  The  station  weight  w-z  is  inversely  proportional  to  the  number  of  traces 
in  30°  azimuthal  bins  to  lessen  the  likelihood  that  waveforms  from  a  cluster  of  closely  spaced 
stations  dominate  the  solution.  The  amplitude  weight  is  inversely  proportional  to  the  time 
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integration  of  the  square  of  the  amplitude  of  observed  waveforms,  so  body-wave  time  windows 
contribute  approximately  equally  to  the  misfit  function. 

Partly  because  waveforms  can  be  shifted  in  time,  the  moment  tensor  solution  and  the 
waveform  misfit  is  insensitive  to  lateral  changes  in  location.  Yet  for  some  traces,  the  onset  times 
of  first  arrivals  can  be  picked  accurately.  To  integrate  this  arrival  time  information  into  the 
global  optimization  of  the  source  solution,  we  calculate  a  transformed  misfit  function  by 
multiplying  the  weight  w,  with  an  additional  first-arrival-time-related  factor  after  the  moment 
inversion.  This  posterior  factor  does  not  affect  the  moment  inversion  at  individual  grid  but  helps 
to  select  the  best  amplitude-fitting  solutions  close  to  the  best  location(s)  required  by  the  travel 
times  of  the  reference  model.  Unlike  other  weighting  factors,  it  is  a  function  of  grid  search 
location  wf(p)  =  (l  +  \Ati(p)\)2/(l  +  aTi)2,  where  At  is  the  time  difference  between  the  onset 

times  of  the  observed  and  synthetic  first  arrivals,  T  the  travel  time  of  the  first  arrival,  and  a  the 
scaling  factor  (-0.01)  to  account  for  a  larger  uncertainty  in  At  at  a  greater  distance  due  to  errors 
in  the  velocity  model.  For  traces  with  no  picked  first  arrival  time,  wf'(p)  =  1.  For  noise-free 
synthetics,  the  picking  of  the  onset  of  the  first  arrival  involves  the  detection  of  the  first  value 
deviated  from  zero  (above  numerical  errors),  which  is  trivial.  The  onsets  of  the  observed  first 
arrivals  are  picked  by  hand  (though  an  automated  first-arrival  picking  scheme  can  be 
implemented  for  routine  and  quick  processing),  with  the  synthetics  generated  from  the  SGTs  and 
an  explosive  source  at  the  EHB  earthquake  location  as  guidance  for  the  pick.  Since  the  SGTs  are 
saved  at  1  s  time  interval,  we  estimate  that  At  has  an  uncertainty  of  -  1  s  (For  a  study  just  on 
source  mechanisms,  one  may  save  the  SGTs  at  much  smaller  time  intervals  at  the  target 
locations).  The  traces  with  first  arrival  times  thus  have,  as  a  group,  relatively  larger  weights.  As 
grid  search  moves  away  from  the  best  travel-time  solution(s),  the  traces  with  first-arrival  times 
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have  even  larger  weights.  This  effectively  panalizes  solutions  that  predict  the  onsets  of  first 
arrivals  significantly  different  from  the  observed  and  causes  the  transformed  misfit  to  increase 
more  quickly  away  from  the  best  solution  than  that  giving  no  consideration  to  the  onset  times  of 
first  arrivals.  For  most  earthquakes  in  this  study,  this  posterior  treatment  of  the  global  misfit 
function  does  not  affect  the  moment  tensors,  but  changes  the  estimates  of  the  location  errors. 

3.4  Inversion  procedures 

For  a  given  source  location,  the  synthetics  are  linear  combinations  of  the  elementary  moment 
tensors  (equation  3).  But  waveform  varies  in  a  nonlinear  fashion  as  the  source  location  changes. 
To  deal  with  this  linear-nonlinear,  moment-tensor-and-location  problem,  we  adopt  a  grid-search 
approach,  solving  the  linear  problem  at  each  grid  in  a  3D  volume  surrounding  the  reference 
(EHB)  source  location.  The  SGT  databases  are  well  suited  for  this  approach,  as  this  involves 
only  extracting  SGTs  at  the  target  3D  grids.  We  limit  the  search  to  a  volume  of  approximately 
91  km  x  91  km  x  28  km  in  the  latitude,  longitude,  and  vertical  directions,  respectively,  centered 
on  the  reference  source  location.  This  translates  to  7  x  7  x  15  grid  points,  or  735  moment  tensor 
solutions,  one  for  each  grid.  Adjustments  are  made  for  events  with  the  reference  locations  near 
the  surface  (no  search  above  the  free  surface).  Because  of  discrete  search,  the  solution  has  an 
intrinsic  numerical  location  error  of  approximately  half  grid  (6.5  km  laterally  and  1  km 
vertically),  which  is  comparable  to  or  smaller  than  the  location  error  caused  by  inaccurate 
reference  models  and  the  error  in  the  earthquake  catalog.  As  the  velocity  model  improves,  it 
may  warrant  saving  the  SGTs  at  every  finite-difference  wave  simulation  grid  (~4.4  km  laterally 
and  1  km  vertically)  or  finer  grid  (by  interpolation)  in  the  target  volumes  of  earthquakes. 

At  each  grid,  we  obtain  an  initial  solution  from  the  normal  equations  using  only  longer- 
period  waves  (30-60  s  and  60-120  s).  Only  the  basis  displacements  p)  for  the  elementary 
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moment  tensors  are  used  (equations  6  and  7)  and  no  synthetic  waveforms  are  required  at  this 
stage.  For  long-period  waves,  a  small  phase  shift  is  less  problematic.  Because  the  earthquakes 
studied  are  all  tectonic  events  as  opposed  to,  say,  explosions,  we  solve  for  pure-deviatoric 
moment  tensors,  so  the  number  of  unknown  coefficients  for  the  elementary  moment  tensors  is  5. 
The  deviatoric  moment  tensor  is  constructed  from  the  coefficients  as  in  Kikuchi  &  Kanamori 
(1991).  After  we  obtain  the  initial  solution,  we  construct  the  synthetics  using  equation  3.  The 
time-windowed  and  bandpass-filtered  synthetics  are  cross-correlated  with  the  observed,  and  the 
synthetics  are  shifted  by  the  time  lag  of  maximum  cross-correlation.  We  obtain  a  new  moment 
tensor  solution  using  the  time-shifted  synthetics.  Since  a  new  moment  tensor  yields  somewhat 
different  synthetics  and  different  time-shifts,  this  process  is  repeated  until  the  solution  is  stable. 
We  find  the  solution  usually  becomes  stable  after  2  iterations. 

After  we  loop  through  all  the  grids,  the  solution  (the  moment  tensor  and  grid  location)  that 
yields  the  minimum  misfit  is  selected  as  our  tentative  preferred  solution.  We  then  check  the 
cross -correlation  coefficients  between  the  bandpass-filtered,  time  windowed  and  shifted 
synthetics  from  the  tentative  solution  and  the  observed.  If  all  the  cross-correlation  coefficients 
are  greater  than  a  prescribed  value  (0.5),  then  the  tentative  solution  becomes  the  final  preferred 
solution.  Otherwise  a  second  and  final  round  of  grid  search  is  carried  out  with  the  frequency¬ 
time  windows  having  cross-correlation  coefficient  less  than  the  cutoff  value  removed  from  the 
inversion.  The  poor  waveform  similarity  for  the  best  solution  in  the  first  round  of  grid  search 
can  be  due  to  (1)  non-stationary  random  noise  that  is  stronger  in  the  time  window  than  before  the 
first  arrival  (SNRs  are  calculated  from  the  noise  level  before  the  first  arrival),  and  (2) 
heterogeneities  that  are  not  accounted  for  in  the  3D  reference  models.  In  either  case,  if  poor- fit 
frequency-time  windows  are  not  removed  from  the  inversion,  they  may  significantly  degrade  the 
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moment  tensor  solution.  Liu  et  al.  (2004)  select  good  data  traces  based  on  synthetic  waveforms 
of  already  good  source  parameters  determined  previously  in  other  studies.  Our  approach 
provides,  in  essence,  the  same  quality  control  on  input  data  without  a  prior  solution. 

We  estimate  the  location  errors  by  examining  the  variation  of  the  misfit  function  within  the 
grids.  We  adopt  the  approach  of  Sambridge  &  Kennett  (1986)  and  normalize  the  misfit  function 
by  the  misfit  of  our  preferred  solution  (the  minimum  misfit).  The  region  that  satisfies  the 
confidence  limit  (1-a)  is 

[E{m)-E(msol)]^Z2Ndfa-a) 

E(  mso1)  Ndf 

where  £(m'"/)  is  the  misfit  of  our  preferred  solution  and  Ndf  the  number  of  degrees  of  freedom 
for  the  chi-squared  distribution.  Since  there  are  five  elementary  moment  tensors  for  deviatoric 
earthquakes  and  three  location  parameters  (latitude,  longitude  and  depth),  the  total  degree  of 
freedom  is  8.  Because  we  allow  the  synthetics  to  shift  in  time  in  waveform  fitting,  we  do  not 
consider  the  origin  time  another  degree  of  freedom,  though  increasing  the  degree  of  freedom  to  9 
does  not  significantly  change  the  right  hand  side.  At  the  one-standard  deviation  (a,  68% 
confidence  limit)  level,  jg(0.68)/8  =  1.159  versus  jg(0.68)/9  =  1.156.  At  the  95%  confidence 
limit,  the  normalized  chi-squared  values  are  1.938  and  1.879,  respectively. 

In  the  following  section  we  present  examples  of  the  source  mechanisms  determined  from  the 
method  outlined  above.  For  events  in  the  GCMT  catalog,  we  compare  the  waveform  fits  for  the 
new  FDSGT  and  GCMT  solutions. 
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4  FDSGT  MOMENT  TENSORS 

The  October  27,  2001  (Ms  5.5)  Yunnan  earthquake  (event  2001.300.05.35.42  in 
Supplementary)  was  recorded  by  the  HIMNT  temporary  network  stations  (de  la  Torre  & 
Sheehan,  2005)  and  the  permanent  stations  at  regional  and  teleseismic  distances  (Fig.  3g).  Body 
and  surface  waves  at  both  regional  and  teleseismic  distances  in  all  frequency  bands  provide 
useful  constraints  on  the  source  parameters,  though  the  SNR  for  this  earthquake  is  best  in  the  30- 
60  s  period  range  (Fig.  3a-d).  Approximately  l/3ld  of  the  frequency-time  windows  are  rejected 
in  the  final  inversion  because  of  their  low  SNRs  or  correlation  coefficients  calculated  before  and 
after  the  first  round  of  grid  search,  respectively. 

The  source  locations  in  the  EHB  and  GCMT  catalogs  differ  by  ~50  km  (Fig.  3f).  For 
moderate  and  small  (Mw<6.5)  earthquakes  like  this  one,  the  difference  in  location  is  unlikely  the 
effects  of  finite  spatial  rupture,  but  dominated  by  location  errors  (Smith  &  Ekstrom,  1997).  Our 
new  solution  is  southeast  of  the  EHB  location  and  southwest  of  the  GCMT  solution,  though  the 
EHB  location  is  within  the  one-a  uncertainty  region  defined  by  waveform  misfit  and  both  the 
EHB  and  GCMT  locations  are  within  the  95%  confidence  limit.  The  depth  of  the  preferred 
solution  (8±3  km)  is  significantly  shallower  (at  the  1-a  uncertainty  level)  than  the  GCMT  depth 
(15  km)  and  the  EHB  depth  (16.7±3.2  km).  The  waveform  fit  at  the  GCMT  and  EHB  depths 
(15-17  km)  is  substantially  worse  than  the  fit  at  the  best  solution  (Fig.  3e).  Compared  to  the 
traveltime-based  EHB  solution  and  long-period  waveform-based  GCMT  solution,  the  depth 
resolution  in  this  study  comes  from  variations  in  waveform  as  a  function  of  the  event  depth  and 
shorter-period  body  waves  used  in  this  study.  The  above  comparison  is  based  on  the  FDSGT 
moment  tensor  solutions  at  the  individual  grids.  If  we  use  the  GCMT  moment  tensor  at  the 
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closest  grid  to  the  GCMT  location,  the  waveform  misfit  is  nearly  three  times  larger  than  that  of 
the  best  FDSGT  solution  (1.64  of  GCMT  versus  0.52  of  FDSGT). 

Figure  4a  shows  the  lowest  misfit  solution  at  each  grid-search  depth.  The  solutions  at  very 
shallow  (2  km)  and  large  (>  18  km)  depths  are  noticeably  different  from  the  global  best  solution. 
In  contrast,  a  ~10  km  lateral  shift  in  location  does  not  affect  the  moment  tensor  solution 
significantly  (Fig.  4b).  This  low  sensitivity  to  lateral  shift  in  location  is  partly  due  to  time-shift 
in  waveform  fitting  to  account  for  errors  in  the  3D  reference  model. 

A  close  comparison  of  the  observed  waveforms  and  the  synthetics  predicted  by  the  GCMT 
and  our  FDSGT  solutions  reveals  that  the  FDSGT  location  yields  waveforms  that  are  better 
matched  in  time.  To  illustrate  this,  we  zoom  in  on  the  first  arrivals  at  station  LSA  and  CHTO, 
and  the  teleseismic  S  wave  at  station  MDJ  (Fig.  5).  At  station  LSA,  the  GCMT  and  FDSGT 
synthetics  essentially  overlap  and  both  are  ~4±1  s  later  than  the  observed.  The  delay  may  result 
from  a  ~3%  lower- than -real  upper  mantle  P  velocity  in  the  reference  model.  Since  a  significant 
portion  of  the  event  to  LSA  path  cuts  through  the  Indian  mantle  indentor,  the  velocity  anomaly 
seems  possible.  At  station  CHTO,  the  onset  of  the  first  arrival  from  the  FDSGT  solution 
matches  fairly  well  with  the  observed.  Changing  the  source  to  the  GCMT  location  (and  the 
GCMT  moment  tensor),  we  have  an  onset  of  the  synthetic  first  arrival  that  appears  ~10  s  later 
than  the  observed.  To  fit  the  arrival,  the  upper  mantle  P  velocity  along  the  event  to  CHTO  path 
has  to  be  8%  higher  than  in  the  reference  model,  a  value  too  large  to  be  likely  true.  At  station 
MDJ,  the  teleseismic  P  wave  is  below  the  SNR  threshold.  The  synthetic  teleseismic  S  wave 
from  the  FDDGT  solution  arrives  slightly  later  than  that  of  the  GCMT  solution  and  fits  the 
observed  better  in  time. 
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Also  shown  in  Figure  5  are  the  synthetics  of  200  random  realizations  of  the  FDSGT  solution. 
In  each  realization,  we  take  each  of  the  estimated  standard  deviations  for  the  six  moment  tensor 
components  and  change  each  component  by  a  random  value  that  follows  a  normal  distribution 
with  the  standard  deviation  for  this  component.  These  realizations  thus  reflect  only  the 
uncertainties  of  the  moment  tensor  components  and  do  not  include  the  uncertainties  of  the  source 
location.  There  are  two  underlying  assumptions  in  this  calculation:  (1)  the  waveform  misfits  are 
due  to  random  noise,  and  (2)  the  moment  tensor  components  are  independent.  These 
assumptions  are  unlikely  to  be  fully  true,  since  part  of  the  waveform  misfit  is  certainly  caused  by 
un-modeled  velocity  heterogeneities  and  thus  not  fully  random.  Nevertheless,  the  range  of 
variations  in  the  200  realizations  is  comparable  to  the  observed  misfit  in  waveforms  after  time 
shift  (also  see  Fig.  3),  suggesting  that  the  estimated  level  of  uncertainties  in  the  moment  tensor 
solution  is  correct.  Increasing  the  moment  tensor  uncertainties  by  a  factor  of  2  results  in 
realizations  that  may  be  opposite  in  phase  polarity  (not  shown).  On  the  other  hand,  reducing  the 
moment  tensor  uncertainties  by  a  factor  of  2  results  in  a  range  of  waveform  variations  that  is 
much  smaller  than  the  observed.  We  note  that  the  middle  of  the  S  wave  train  at  station  MDJ  is 
particularly  sensitive  to  errors  in  the  moment  tensor  compared  to  the  rest  of  the  waveform.  The 
variation  in  synthetic  waveforms  due  to  uncertainties  in  the  moment  tensor  together  with  the 
estimates  of  the  source  location  errors  can  be  used  to  provide  proper  weighting  of  the  data  in 
waveform-based  tomographic  inversion. 

The  June  29,  2002  (Ms  5.5)  Qinghai  earthquake  (event  2002.180.06.54.42  in  Supplementary) 
was  also  recorded  by  the  HIMNT  temporary  network  stations  (de  la  Torre  &  Sheehan,  2005)  and 
the  permanent  stations  at  regional  and  teleseismic  distances  (Fig.  6g).  Our  new  solution  is  ~25 
km  south  of  the  EHB  and  GCMT  locations,  which  are  outside  of  the  95%  confidence  limit  (Fig. 
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6f).  Inspection  of  the  first  arrivals  at  KMI,  LSA  and  the  HIMNT  stations  confirms  that  the 
southward  shift  in  location  from  the  reference  is  required  by  the  data  (and  the  3D  reference 
model).  At  8±3  km  depth,  the  new  solution  is  also  significantly  shallower  than  the  EHB  (16.4 
km)  and  GCMT  (33  km)  estimates  (Fig.  6e).  Despite  the  ~25  km  offset  in  the  lateral  locations, 
the  FDSGT  and  GCMT  moment  tensors  are  very  similar  (Fig.  7),  reflecting  again  the  low 
sensitivity  of  the  moment  tensor  solutions  to  lateral  shift  in  source  location. 

Figure  7  shows  a  comparison  of  the  FDSGT  and  GCMT  solutions  for  13  earthquakes.  For 
most  of  the  earthquakes,  the  two  solutions  are  in  general  agreement.  One  notable  exception  is 
the  August  18,  2003  earthquake  at  the  eastern  Himalayan  syntaxis.  The  GCMT  solution  shows  a 
strike-slip  event  at  33  km  depth,  while  our  FDSGT  solution  is  a  shallow  (2+16  km)  thrust  event. 

No  all  FDSGT  solutions  are  significantly  shallower  than  the  GCMT  solutions.  Four  out  of 
thirteen  are  significantly  shallower  than  the  GCMT  solutions,  one  significantly  deeper,  and  the 
rest  within  the  95%  confidence  limit.  The  one  with  a  FDSGT  source  depth  significantly  deeper 
than  the  GCMT  solution  (the  July  21,  2003  earthquake,  event  2003.202.15.16.37  in 
Supplementary)  has  a  local  minimum  near  the  GCMT  depth  (15  km),  but  the  global  minimum  is 
at  46.7±3.6  km  depth.  This  example  highlights  the  benefit  of  a  global  optimization  approach  in 
earthquake  source  parameter  study,  such  as  grid  search.  Though  if  the  grid  search  volume  is  not 
large  enough  for  a  poor  reference  location,  one  may  still  miss  the  real  solution.  A  perturbation 
approach  about  the  reference  runs  the  risk  that  the  solution  is  trapped  at  a  local  minimum. 

Also  shown  in  Figure  7  are  the  FDSGT  solutions  for  10  additional  earthquakes  with  no 
corresponding  GCMT  solutions,  de  la  Torre  et  al.  (2007)  analyzed  the  December  2,  2001 
earthquake  (2001.336.22.41.14  in  Supplementary)  along  the  Himalayan  arc  in  northeast  India 
and  identified  it  as  a  strike-slip  event  with  P  axis  trending  to  the  south,  a  result  generally 
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consistent  with  our  FDSGT  determination.  A  list  of  the  moment  tensor  solutions  and  additional 
waveform  fits  are  available  in  Supplementary  Materials. 

5  DISCUSSION  AND  CONCLUSION 

We  determined  source  parameters  for  23  earthquakes  between  2001  and  2004  in  the 
southeast  Tibetan  plateau,  using  the  SGTs  of  49  temporary  and  permanent  seismic  stations  at 
regional  and  teleseismic  distances.  A  3D  reference  model  based  on  CUB  2.0  is  used  in  the 
calculation  of  the  SGTs.  Taking  advantage  of  the  SGT  database,  we  search  the  best  source 
location  and  the  corresponding  moment  tensors  in  a  volume  surrounding  the  reference  EHB 
location.  At  each  location,  the  moment  tensor  is  determined  by  minimizing  the  waveform  misfit 
function  of  time-windowed  first  arrivals,  teleseismic  S  waves  and  surface  waves.  The 
waveforms  are  bandpass-filtered  in  four  frequency  bands  to  accommodate  variations  in  the 
quality  of  waveforms  at  different  frequencies.  The  waveforms  are  time  shifted  to  accommodate 
errors  in  the  SGTs  due  to  the  imperfect  reference  model.  No  prior  source  mechanism  solution  is 
needed  in  the  inversion. 

Our  FDSGT  moment  tensor  determinations  are  in  general  agreement  with  the  GCMT 
solutions,  though  there  are  significant  differences  for  some  earthquakes.  The  differences  are 
likely  due  to  differences  in  the  data  used  in  the  GCMT  inversion  and  in  this  study,  in  addition  to 
the  differences  in  the  Green  functions  (ID  versus  3D  reference  models).  Because  waveforms  are 
time-shifted  before  inversion,  the  moment  tensor  solutions  are  not  sensitive  to  uncertainties  in 
lateral  source  location,  and  the  FDSGT  and  GCMT  source  locations  may  differ  up  to  10s  km. 
Inspection  of  the  observed  and  synthetics  in  details  shows  that  regional  stations,  the  3D  reference 
model,  and  shorter-period  body  waves  in  this  study  tend  to  provide  tighter  constraints  on  the 
lateral  position  as  well  as  depth. 
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The  benefits  of  FDSGT  moment  tensor  inversion  are  likely  greatest  for  small  to  moderate 
earthquakes  (e.g.,  the  additional  earthquakes  with  no  corresponding  GCMT  solutions  in  Figure 
7),  which  usually  do  not  generate  long-period  waves  at  teleseismic  distances,  the  records  used  in 
GCMT  solutions.  The  lowest  magnitude  of  the  earthquakes  in  this  study  is  Mw  4.35,  though  the 
method  is  applicable  to  smaller  earthquakes.  We  note  that  the  3D  SGT  approach  is  most 
effective  when  the  model  used  to  represent  earth  structure  predicts  significantly  better  waveform 
fit  than  ID  models.  The  3D  reference  model  used  in  this  study  does  not  adequately  predict  short- 
period  (<15  s)  surface  waves,  so  only  long -period  surface  waves  and  intermediate-period  body 
waves  are  used  in  the  inversion.  As  the  crustal  velocity  model  improves,  one  may  extend  the 
useful  surface  waves  to  shorter  periods  and  be  able  to  determine  the  moment  tensors  of  lower 
magnitude  earthquakes. 
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Figure  1.  Distributions  of  the  broadband  seismic  stations  used  in  source  inversion  (triangles)  and 
the  earthquakes  studied  (dots).  White  rectangles  outline  the  wave  simulation  models  at  (1) 
regional,  (2)  far-regional,  and  (3)  teleseismic  scales. 
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Figure  2.  Phase  delays  and  maximum  cross-correlation  coefficients  between  the  observed  and 
synthetic  waveforms  at  four  frequencies.  The  CRUST  2.0  model  predicts  a  larger  average  phase 
delay  and  lower  maximum  cross -correlation  coefficient  (dashed  lines)  than  the  CUB  2.0  model 
(solid  lines).  The  individual  CRUST  2.0  values  (circles)  are  also  more  scattered  than  the  CUB 
2.0  predictions  (crosses).  The  GCMT  moment  tensor  for  this  August  24,  2004  earthquake  is 
used  in  the  calculation  of  the  synthetics. 
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Figure  3.  (a-d)  Bandpass -filtered  observed  waveforms  (thick  gray  lines)  are  compared  with  the 
synthetics  for  the  best  source  solution  (red  lines).  The  frequency  bands  are,  from  left  to  right, 
0.00833-0.0167,  0.0167-0.0333,  0.0333-0.0667,  and  0.0667-0.167  Hz.  Dark  thick  lines  are  time- 
shifted  observed  waveforms  used  in  the  final  inversion.  For  the  first  arrivals  (in  c  and  d),  the 
time  is  reduced  to  align  the  arrivals  within  the  time  window.  The  amplitude  is  normalized  by 
individual  trace,  (e)  Normalized  waveform  misfit  as  a  function  of  depth.  Crosses  mark  the  grid- 
search  depths.  The  shaded  region  indicates  the  1-a  confidence  limit,  (f)  Contours  of  various 
confidence  limits  at  the  depth  of  the  global  best  solution  (8  km  for  this  event).  The  cross  marks 
the  FDSGT  location,  the  square  the  EHB  location,  and  the  diamond  the  GCMT  location.  From 
the  innermost  contour  surrounding  the  FDSGT  solution,  the  contours  are  the  68%  (1-a,  dark 
line),  80%,  90%,  95%  (2-a,  dark  line),  and  99%  confidence  limits,  (g)  Triangles  and  cross  mark 
the  stations  used  in  the  source  inversion  and  the  event,  respectively.  Gray  contours  follow  the 
sea  level,  and  2000  and  4000  m  elevation. 
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Figure  4.  Variations  of  the  best  moment  tensor  solution  as  a  function  of  depth  (upper  panel)  and 
lateral  locations  at  the  global  best  solution  depth  (lower  panel).  The  global  best  solution  is  shown 
as  the  blue-and-white  beach  ball.  The  GCMT  solution  (red- and- white  beach  ball)  is  drawn  for 
comparison.  The  background  is  shaded  topographic  relief  with  known  fault  traces  marked  by 
dashed  lines. 
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Figure  5.  A  comparison  of  the  observed  waveforms  (thick  gray  lines)  and  the  synthetics  for  the 
FDSGT  (solid)  and  GCMT  (dashed)  solutions  at  stations  LSA,  CFITO  and  MDJ  (triangles).  The 
200  random  realizations  of  the  FDSGT  solution  are  plotted  below  the  observed  waveforms  to 
show  the  range  of  waveform  variation  due  to  uncertainties  in  the  moment  tensor  solution.  The 
waveforms  have  been  filtered  between  0.01667-0.0333  Hz,  but  not  shifted  in  time.  The  cross 
and  circle  near  26°  latitude  and  101°  longitude  are  the  FDSGT  and  GCMT  locations, 
respectively.  The  channel  and  time  window  (after  the  earthquake  origin)  are  marked  above  the 
waveforms. 
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Figure  6.  Same  as  Fig.  3,  except  for  the  June  29,  2002  earthquake 
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Figure  7.  A  comparison  of  the  FDSGT  (white-and-blue  symbols)  and  GCMT  (white-and-red 
symbols)  solutions  for  13  earthquakes  in  the  southeast  Tibetan  plateau.  The  numbers  above  the 
FDSGT  solution  denote  the  year  and  Julian  day  of  the  earthquake.  Also  plotted  are  the  source 
mechanisms  of  10  additional  earthquakes  without  GCMT  solutions  (black-and-white  symbols). 
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SUPPLMENTARY  MATARIALS 

Moment  tensor  solutions  based  on  finite-difference  stain  Green  tensor  databases 
The  following  list  contains  the  source  locations  and  moment  tensors  of  23  earthquakes  in  the 
southeast  Tibetan  plateau.  It  uses  two  lines  for  each  earthquake. 


First  line:  Event  ID,  location,  time  and  moment  magnitude.  Numbers  stand  for  the  positions  of 
the  variables  in  the  line,  which  are  separated  by  space. 

I  Event  ID  in  the  format  of  year.J  ulian-day.hour.minite.second 
2-4  EHB  latitude  and  longitude  in  degree,  and  depth  in  kilometer 
5-6  Latitude  and  one  standard  deviation  in  degree 

7-8  Longitude  and  one  standard  deviation  in  degree 
9-10  Depth  and  one  standard  deviation  in  kilometer 

II  Time  with  respect  to  the  reference  time  in  second 
12  M  w  as  defined  in  Dahlen  and  T romp  (1998,  eq.5.91) 


Second  line:  Moment  tensor  components 

1-12  The  six  moment-tensor  elements:  M  rr,  M  tt,  M  pp,  M  rt,  M  rp,  M  tp,  each  followed  by  its  one 
standard  error.  The  subscript  r  denotes  up,  t  south,  and  p  east. 

13  The  exponent  shared  by  all  the  moment-tensor  elements  and  their  standard  errors.  An 
exponent  of  23  stands  for  1023  dyne-cm. 

2001.300.05.35.42  26.295  100.578  16.7  26.240  0.120  100.710  0.068  8.00  3.05  -1.38  5.39 

-1.81e+00  1.54e+00  -1.37e+01  2.75e+00  1.55e+01  2.28e+00  2.38e+00  1.07e+00  8.01e-01 1.57e+00  -4.37e+00  1.99e+00  23 
2001.336.22.41.14  27.281  88.236  20.9  27.440  0.000  88.290  0.068  6.00  6.55  -3.18  4.83 

1.07e+00  2.02e-01  -2.56e+00  2.82e-01 1.49e+00  1.97e-01  -1.79e-01  3.91e-01 1.43e-01 1.64e-01  -1.30e-01 1.15e-01  23 
2002.180.06.54.42  34.043  94.393  16.4  33.800  0.060  94.500  0.270  8.00  4.05  -1.11  5.45 

-7.15e+00  8.55e-01  4.93e+00  1.00e+00  2.22e+00  5.19e-01  -8.56e-01  8.81e-01  -1.37e+00  3.84e-01 1.75e+01  8.57e-01  23 
2002.220.11.42.07  30.862  99.889  35.0  30.560  0.000  100.305  0.068  16.64  1.23  -0.21  4.99 

-3.09e+00  7.71e-02  4.17e+00  1.15e-01  -1.08e+00  8.51e-02  -4.55e-01  7.45e-02  9.03e-03  2.73e-02  4.83e-01  3.45e-02  23 
2002.270.17.14.40  33.324  93.484  35.0  32.960  0.000  93.690  0.068  50.34  22.76  1.82  4.93 

-1.04e+00  2.33e-01  -4.69e-01  4.67e-01 1.51e+00  4.04e-01  -1.38e+00  2.61e-01  -1.24e-01 1.01e-01  2.43e+00  4.77e-01  23 
2003.202.15.16.37  25.990  101.288  35.0  25.760  0.000  101.115  0.067  54.14  7.22  -0.74  5.91 

4.09e+01  4.00e+01  -8.24e+01  5.29e+01  4.15e+01  3.46e+01  4.05e+01  4.43e+01  -1.30e+01  2.09e+01  -4.23e+01  3.12e+01  23 
2003.220.09.59.24  26.738  97.003  6.4  26.960  0.000  97.200  0.202  6.00  5.05  0.26  4.70 

7.32e-01  6.30e-02  -3.21e-01  6.76e-02  -4.11e-01  2.47e-02  -1.05e+00  1.24e-01  4.97e-01  7.33e-02  4.50e-01  2.95e-02  23 
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2003.230.09.03.02  29.564  95.574  16.1  29.480  0.120  95.850  0.338  2.00  15.60  -5.34  5.49 

6.41e+00  1.32e+00  -9.35e+00  1.41e+00  2.94e+00  4.90e-01  -1.98e+01  2.55e+01  3.74e+00  3.35e+00  -1.09e+00  7.60e-01  23 
2003.238.20.54.17  24.973  93.619  45.3  25.160  0.120  93.555  0.337  46.67  23.60  0.34  4.92 

6.51e-01  2.14e-01  -2.78e+00  3.47e-01  2.13e+00  2.74e-01  -1.26e-01  9.51e-02  6.26e-01  7.26e-02  1.66e+00  2.73e-01  23 
2003.289.12.28.11  25.917  101.234  35.0  25.760  0.060  101.250  0.135  50.34  11.23  -2.90  5.56 

-5.03e+00  2.13e+00  -1.80e+01  4.21e+00  2.30e+01  3.63e+00  3.13e+00  9.80e-01  -4.27e+00  9.48e-01  -1.68e+01  2.51e+00  23 
2003.330.13.39.00  27.311  103.726  35.0  27.440  0.000  103.410  0.068  54.14  5.91  0.07  4.87 

5.20e-01 1.95e-01  9.96e-01  3.27e-01  -1.52e+00  2.63e-01  5.68e-01 1.25e-01  -8.00e-01 1.47e-01  -1.88e+00  2.99e-01  23 

2003.336.23.26.51  30.816  101.078  35.0  31.160  0.060  101.385  0.135  62.16  17.48  -1.73  5.01 

1.44e-01  2.31e-01  2.83e+00  3.53e-01  -2.98e+00  2.67e-01  7.81e-01 1.57e-01 1.35e+00  1.09e-01  -2.41e+00  4.94e-01  23 

2003.364.20.59.49  25.381  96.469  35.0  25.640  0.060  96.660  0.068  19.09  13.24  -1.14  4.35 

-6.24e-02  5.20e-02  3.34e-01 1.02e-01  -2.71e-01  8.76e-02  1.47e-01  2.77e-02  2.19e-03  1.38e-02  2.39e-01  8.17e-02  23 

2004.050.09.11.42  25.473  96.776  35.0  25.760  0.060  96.390  0.338  46.67  22.76  -2.03  4.73 

-2.25e-01  7.07e-02  1.59e+00  2.65e-01  -1.37e+00  2.56e-01  2.97e-01  5.67e-02  -8.84e-02  2.56e-02  2.58e-01  3.13e-01  23 
2004.066.11.54.43  33.298  91.892  35.0  33.080  0.000  92.070  0.203  30.13  6.01  -0.66  4.85 

5.43e-01  2.90e-01  -1.61e-01  5.38e-01  -3.82e-01  4.54e-01  -1.35e+00  1.74e-01  -1.08e+00  2.78e-01 1.59e+00  3.30e-01  23 
2004.067.13.29.44  31.581  91.153  4.4  31.640  0.000  91.395  0.068  10.00  3.05  -0.48  5.30 

-4.02e+00  2.03e+00  -7.24e+00  2.61e+00  1.13e+01 1.64e+00  1.47e-01 1.88e+00  -3.87e+00  1.43e+00  3.49e+00  8.68e-01  23 
2004.085.21.15.01  25.953  101.220  35.0  25.760  0.060  101.385  0.067  50.34  11.23  -0.61  4.82 
-5.55e-01 1.51e-01  -8.91e-03  3.27e-01  5.64e-01  2.90e-01 1.24e-01 1.04e-01  3.57e-01  7.37e-02  -2.04e+00  3.68e-01  23 
2004.087.18.47.29  33.912  89.128  10.7  33.800  0.060  89.370  0.270  8.00  2.00  1.48  5.83 

-6.05e+01  2.02e+01  -4.49e-01  2.54e+01  6.09e+01 1.54e+01 1.31e+01  2.47e+01  -3.02e+01 1.99e+01 1.15e+01  6.64e+00  23 
2004.171.10.24.53  25.264  95.499  35.0  25.520  0.180  95.175  0.202  24.37  18.75  0.24  4.41 

1.32e-01  2.22e-02  -1.95e-01  3.34e-02  6.37e-02  2.49e-02  -2.46e-01 1.28e-02  -1.55e-01  7.80e-03  -3.92e-01  3.05e-02  23 
2004.185.14.10.48  34.054  89.276  35.0  33.680  0.060  89.370  0.202  16.64  22.76  -1.63  4.99 

-2.41e+00  2.06e-01  -1.36e+00  3.76e-01  3.77e+00  3.15e-01  3.54e-01 1.83e-01  -1.80e+00  4.05e-01  9.88e-01  2.17e-01  23 
2004.205.01.25.31  30.080  88.010  4.9  30.440  0.120  88.290  0.068  8.00  3.00  -0.94  4.50 

-3.87e-01  3.82e-01  -1.17e-01  4.18e-01  5.04e-01 1.69e-01  4.70e-01  5.72e-01  -2.62e-01 1.52e-01 1.47e-01  2.27e-01  23 
2004.223.10.26.15  27.310  103.887  12.6  27.320  0.000  103.545  0.067  19.09  9.02  1.71  5.19 

2.85e+00  1.62e+00  1.07e+00  2.16e+00  -3.92e+00  1.42e+00  -1.53e+00  6.18e-01  -5.65e+00  8.81e-01  -3.71e+00  8.63e-01  23 
2004.237.10.05.33  32.537  92.137  3.4  32.480  0.060  92.205  0.135  8.00  5.32  -1.18  5.30 

-8.97e+00  8.16e-01  -1.98e+00  1.23e+00  1.10e+01  9.17e-01 1.67e+00  7.77e-01  4.82e-01 1.28e+00  4.83e+00  3.82e-01  23 
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Figure  SI.  Same  as  Fig.  3,  except  for  the  J  uly  21,  2003  earthquake  (2003.202.15.16.37). 
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Figure  S2.  Same  as  Fig.  3,  except  for  the  October  16,  2003  earthquake  (2003.289.12.28.11). 
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Figure  S3.  Same  as  Fig.  3,  except  for  the  M  arch  7,  2004  earthquake  (2004.067.13.29.44). 
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Figure  S4.  Same  as  Fig.  3,  except  for  the  M  arch  28,  2004  earthquake  (2004.087.18.47.29). 
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Figure  S5.  Same  as  Fig.  3,  except  for  the  August  24,  2004  earthquake  (2004.237.10.05.33). 
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SUMMARY 

Seismic  tomography  has  been  one  of  the  primary  tools  to  image  the  interior  of  the  earth 
and  other  elastic  structures.  To  date  the  inversions  of  compressional  (P)  and  shear  ( S )  wave 
speeds  have  been  carried  out  separately  under  the  assumption  that  P  traveltimes  are  affected 
only  by  the  P  wave  speed  of  the  elastic  media  and  S  traveltimes  by  the  S  wave  speed.  Using 
numerical  and  analytical  solutions,  we  show  that  for  finite-frequency  seismic  waves,  S  wave 
speed  perturbations  may  have  significant  effects  on  P  waveforms.  This  suggests  that  when 
waveform-derived  traveltime  and  amplitude  anomalies  are  used  in  tomographic  inversions,  the 
P- wave  measurements  should  be  related  to  not  only  P  wave  speed  perturbations  but  also  S 
wave  speed  perturbations. 

Key  words:  Seismic  tomography;  Computational  seismology;  Theoretical  seismology;  Wave 
scattering  and  diffraction;  Wave  propagation. 


INTRODUCTION 

For  many  years,  seismologists  have  made  fundamental  discoveries 
about  the  earth’s  interior,  using  rays  to  describe  the  propagation  of  P 
and  S  seismic  waves.  However,  with  the  increasing  need  for  more  ac¬ 
curate  images  of  the  earth  structure  at  fine  resolution,  the  limitations 
of  ray  theory  become  problematic.  Several  recent  studies  have  de¬ 
veloped  a  new  theory  that  describes  the  sensitivities,  or  Frechet  ker¬ 
nels,  of  the  traveltimes  and  amplitudes  of  finite-frequency  seismic 
waves  to  perturbations  in  the  properties  of  wave-propagation  media 
(Dahlen  et  al.  2000;  Zhao  et al.  2000,  2005;  Tromp  et  al.  2005).  For 
a  finite-frequency  wave,  the  sensitivities  to  the  elastic  and  anelastic 
properties  are  in  general  distributed  in  a  3D  volume  surrounding  the 
ray  path.  Numerical  experiments  have  demonstrated  that,  within  the 
limit  of  Born  approximation,  the  finite-frequency  theory  represents 
the  propagation  of  realistic  seismic  waves  more  accurately  than  ray 
theory  when  the  scales  of  wave  speed  heterogeneities  are  smaller 
than  the  Fresnel  zones  of  the  waves  (Hung  et  al.  2000;  Baig  et  al. 
2003;  Yang  &  Hung  2005;  Zhang  et  al.  2007). 

To  date  the  inversions  of  P  and  S  wave  speeds  have  been  carried 
out  separately  under  the  assumption  that  P  traveltimes  are  affected 
only  by  the  P  wave  speed  of  the  elastic  media  and  S  traveltimes  by 
the  S  wave  speed  (e.g.  Dziewonski  &  Anderson  1984;  Grand  1987; 
Zhao  et  al.  1992;  van  der  Hilst  et  al.  1997).  With  the  approxima¬ 
tions  of  only  forward  scattering  of  body  wave  and  far  held,  Dahlen 
et  al.  (2000)  considered  that  only  P-to-P  scattering  off  a  compres¬ 
sional  wave  speed  heterogeneity  and  like-type  S-to-S  scattering  off 
a  shear  wave  speed  heterogeneity  are  significant.  Although  several 
studies  have  documented  the  scattering  of  wave  held  by  an  arbitrary 
heterogeneity  (Wu  &  Aki  1985)  and  have  calculated  the  sensitivity 
kernels  with  the  near-held  terms  or  in  a  full-wave  approach  (Favier 
et  al.  2004;  Tromp  et  al.  2005;  Zhao  et  al.  2005;  Liu  &  Tromp 


2006) ,  the  assumption  that  the  propagation  of  P  or  5  waves  depends 
solely  on  the  P  or  5  wave  speed  structure,  respectively,  has  remained 
unchanged  in  the  applications  of  the  hnite-frequency  theory  to  date 
(Hung  et  al.  2004;  Montelli  etal.  2006;  Yang  et  al.  2006;  Chen  et  al. 

2007) .  In  this  paper,  we  document  the  contributions  of  S  wave  speed 
perturbations  to  the  P-wave  traveltimes  and  amplitude  anomalies 
using  numerical  and  analytical  solutions.  We  conclude  that  neglect¬ 
ing  the  cross-dependence  between  P  waveforms  and  S  wave  speed 
perturbations  distorts  the  forward  relation  between  the  wave  prop¬ 
agation  medium  and  observations  and  may  cause  a  systematic  bias 
in  the  solution  of  P  wave  speeds. 

MODELS  AND  RESULTS 


To  demonstrate  the  cross-dependence  of  P  and  S  waveforms  and 
wave  speed  perturbations,  we  carry  out  numerical  experiments  us¬ 
ing  a  3-D  staggered-grid,  finite-difference  method  (Olsen  1994). 
We  use  a  simple  homogenous  model  as  our  reference  model,  so  the 
results  can  be  compared  to  and  also  analysed  with  the  analytical 
solutions.  The  P  and  5  wave  speeds  in  the  reference  model  are  6.5 
and  3.5  km  s"1,  respectively.  The  grid  spacing  is  200  m,  and  the 
time  step  0.01  s.  We  introduce  a  cylindrical-shaped  velocity  pertur¬ 
bation  with  a  radius  of  3  km,  a  height  of  6  km  at  various  depths 
along  a  vertical  line  between  the  source  and  receiver,  the  axis  of 
which  is  always  oriented  vertically  in  the  Z  direction  (Fig.  1 ).  The 
traveltime  and  amplitude  anomalies  caused  by  the  wave  speed  per¬ 
turbation  are  measured  by  cross-correlation  of  the  waveforms  with 
and  without  the  wave  speed  perturbation.  For  simplicity,  an  explo¬ 
sive  source  with  a  dominant  period  of  1.2  s  is  used  to  observe  the 
direct  P  arrival.  When  only  the  P  or  S  wave  speed  is  perturbed,  the 
sensitivities  to  the  P  or  S  wave  speed  perturbations  (Fig.  2)  cal¬ 
culated  using  the  scattering-integral  method  (Appendix  A)  predict 
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Figure  1.  The  geometry  of  the  homogeneous  and  isotropic  model  used  in 
the  numerical  simulations.  The  rectangle  with  solid  lines  shows  the  vertical 
cross-section  of  the  model  containing  the  source  (star)  and  receiver  (trian¬ 
gle).  The  horizontal  distance  from  the  source  to  the  receiver  is  53  km.  or 
approximately  7.4  wavelengths  of  the  dominant  P  wave.  The  vertical  line 
in  between  the  source  and  receiver  shows  where  a  cylindrical-shaped  wave 
speed  anomaly  is  placed. 


the  P  traveltime  and  amplitude  anomalies  that  accurately  match  the 
direct  waveform  cross-correlation  measurements  (Fig.  3).  Since  the 
density  remains  the  same  in  all  the  calculations,  a  perturbation  in 
only  the  P  wave  speed  is  equivalent  to  a  perturbation  in  the  bulk 
modulus,  while  a  perturbation  in  only  the  S  wave  speed  corresponds 
to  changes  in  both  the  bulk  and  shear  moduli.  When  both  the  P 
and  S  wave  speeds  are  perturbed  within  the  cylindrical  anomaly, 
however,  the  predicted  P  traveltime  and  amplitude  anomalies  from 
the  sensitivities  to  P  wave  speed  perturbations  alone  do  not  match 
the  direct  waveform  cross-correlation  measurements  (Fig.  3).  Only 
after  we  add  the  contributions  from  the  sensitivities  of  the  P  arrival 
to  5  wave  speed  perturbations,  do  the  kernel  predictions  match  the 
values  from  the  direct  waveform  cross-correlation,  indicating  that 
5  wave  speed  perturbations  affect  the  P  arrival.  On  the  other  hand, 
P  wave  speed  perturbations  have  little  effect  on  the  direct  S  arrival 
in  the  experiment  using  a  double-couple  source  to  generate  S  waves 
(Fig.  3). 

The  sensitivities  of  P  traveltime  to  P  wave  speed  perturbations 
,  Fig.  2)  are  similar  to  the  classic  ‘banana — doughnut'  kernels 
of  Dahlen  et  al.  (2000),  with  a  minimum  sensitivity  along  the  ray 
path,  though  in  this  study  the  ray  path  is  a  straight  line  and  the  kernels 
are  calculated  from  the  full  wavefield.  In  contrast,  the  sensitivities 
of  P  traveltime  to  S  wave  speed  perturbations  are  at  the  maximum 
along  the  ray  path  (Kp  p).  The  patterns  of  the  sensitivities  for  P 
amplitude  variation  to  P  and  S  wave  speed  perturbations  are  also 
reversed:  the  maximum  amplitude  sensitivities  along  the  ray  path 
for  P  wave  speed  perturbations  (K*  )  and  minimum  sensitivities 
along  the  ray  path  for  S  wave  speed  perturbations  (Kpq). 

In  a  homogeneous  and  isotropic  medium,  as  in  the  above  nu¬ 
merical  experiments,  the  sensitivity  kernels  to  P  and  S  wave  speed 
perturbations  can  also  be  calculated  from  the  analytical  wavefield 
solutions  (Aki  &  Richards  2002)  and  the  results  confirm  the  finite- 
difference  calculations.  Using  various  combinations  of  the  near- 
and  far-held  terms  for  the  source  and  receiver,  we  find  that  both  the 
far-  and  near-Held  terms  from  the  source  and  receive  contribute  to 
the  P  traveltime  and  amplitude  sensitivities  to  5  wave  speed  per¬ 
turbations.  In  particular,  the  near-  and  mid-Held  terms  are  mainly 
responsible  for  the  sensitivities  along  the  ray  path  (more  in  the  fol¬ 
lowing  section).  Furthermore,  the  ratio  of  the  maximum,  absolute 


traveltime  sensitivities  to  S  wave  speed  and  P  wave  speed  perturba¬ 
tions  (K%  /K%  )  is  approximately  4kPf)2/(la2)  at  the  mid  point 
between  the  source  and  receiver,  where  kP  is  the  wavelength  of  the 
P  wave,  /  the  source-receiver  distance,  a  the  P  wave  speed,  and  f) 
the  S  wave  speed.  Using  the  parameters  in  our  numerical  experi¬ 
ments  (Fig.  1),  we  find  a  ratio  of  about  0.17,  consistent  with  the 
sensitivity  kernels  (Fig.  2)  and  direct  waveform  cross-correlation 
measurements  (Fig.  3). 

Although  a  homogeneous  and  isotropic  medium  is  an  oversim¬ 
plification  of  the  real  Earth,  the  above  formula  for  the  relative  mag¬ 
nitude  of  the  sensitivities  to  P  and  S  wave  speed  perturbations  pro¬ 
vides  a  simple  way  to  assess,  to  the  first  order,  the  condition  under 
which  the  contributions  of  S  wave  speed  perturbations  to  P  wave¬ 
forms  are  nontrivial.  For  a  medium  with  a  a/fi  of  1. 7-2.0,  a  value 
suitable  for  most  of  the  solid  earth,  the  ratio  of  the  maximum  travel¬ 
time  sensitivities  to  S  wave  speed  and  P  wave  speed  perturbations  at 
the  mid  point  between  the  source  and  receiver  decreases  from  ~0. 17 
for  a  source-receiver  geometry  with  l/kP  of  ~7.4  to  ~0.017  for  a  P 
wave  of  the  same  wavelength  at  10  times  the  propagation  distance 
(l/Xp  =  74). 

To  assess  the  global  relative  contributions  of  P  and  S  wave  speed 
perturbations  to  P  traveltimes  and  amplitudes,  we  integrate  the  ab¬ 
solute  values  of  the  sensitivities  over  the  entire  volume.  For  P- wave 
traveltimes,  the  ratio  of  the  integrated  absolute  values  (Kp  / Kpa  ) 
varies  from  4. 1  per  cent  at  l/k P  of  ~600  to  3 1  per  cent  at  l/'k Pof~l 
(Fig.  4).  For  P-  wave  amplitude  anomalies,  the  ratio  of  the  integrated 
absolute  values  (Kp/Kp)  changes  from  2.7  to  49  per  cent  at  the 
same  corresponding  distances. 


DISCUSSION 

In  the  above  calculations,  the  magnitude  of  5  wave  speed  perturba¬ 
tions  is  the  same  as  the  magnitude  of  P  wave  speed  perturbations. 
In  the  real  Earth,  S  wave  speed  variations  are  usually  larger  than 
P  wave  speed  variations  in  percentage  (e.g.  Karato  &  Karki  2001). 
One  per  cent  of  partial  melt  in  the  upper  mantle,  for  example,  is  esti¬ 
mated  to  cause  3.6  per  cent  P  wave  speed  reduction  and  7.9  per  cent 
5  wave  speed  reduction  (Hammond  &  Humphreys  2000).  The  net 
effect  of  a  larger  5  wave  speed  perturbation  than  the  P  wave  speed 
change  is  an  increase  of  the  contribution  and  importance  of  5  wave 
speed  variations  to  P  waveforms.  P  and  5  wave  speed  variations 
are  also  often  coupled  in  the  real  earth.  In  such  a  case,  a  small- 
scale  perturbation  along  the  ray  path  affects  P  waveforms  when 
both  P  and  S  wave  speeds  are  taken  into  consideration.  This  proves 
to  be  important  in  the  numerical  validation  of  the  sensitivity  kernels 
(Fig.  3). 

Separating  the  far-field  term  (with  l/r  decay,  where  r  is  the  dis¬ 
tance  to  the  source  or  receiver)  from  the  mid-  and  near-field  terms 
(with  l/r2  decay  and  l/r3  decay,  respectively)  in  the  analytical  so¬ 
lution  of  a  homogeneous  and  isotropic  medium  (Aki  &  Richards 
2002),  we  find  that  the  sensitivities  of  P  traveltime  to  P  wave  speed 
perturbations  calculated  from  the  far-field  alone  are  the  same  as 
the  ‘banana-doughnut’  kernels  of  Dahlen  et  al.  (2000),  with  a  zero 
traveltime  sensitivity  along  the  ray  path  (Fig.  5).  The  far-field  P  trav¬ 
eltime  sensitivity  to  5  wave  speed  perturbations  is  also  zero  along 
the  ray  path.  However,  off  the  ray  path  the  sensitivity  is  nontrivial 
compared  to  the  sensitivity  to  P  wave  speed  perturbations  (Fig.  5). 
Unlike  the  scattering  of  a  P  wave  by  a  P  wave  speed  perturbation, 
the  far-field  Rayleigh  scattering  coefficient  of  a  scattered  P  wave 
due  to  an  S  wave  speed  perturbation  is  zero  on  the  ray  path  (fig.  4 
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Figure  2.  (a)  The  sensitivity  kernels  for  the  horizontal  component  of  the  direct  P  arrival  in  the  vertical  plane  containing  the  source  and  receiver.  and  K% 
are  the  traveltime  and  amplitude  sensitivity  kernels  to  P  wave  speed  perturbations,  respectively;  Kp  and  Kp  are  the  traveltime  and  amplitude  sensitivities 
to  S  wave  speed  perturbations,  respectively.  The  unit  of  the  traveltime  sensitivity  is  s  m~3  and  that  of  the  amplitude  sensitivity  is  1  m-3  (Zhang  et  al.  2007). 
(b)  The  K% p,  q,  Kpp  and  Ksp  sensitivities  at  three  vertical  planes  perpendicular  to  the  ray  path.  The  black  line  crossing  the  three  planes  represents  the 
ray  path.  The  source  is  located  at  the  left-hand  side  of  this  figure.  The  three  planes  are  6.5,  26.5  and  46.5  km  away  from  the  source,  respectively.  Less  saturated 
and  different  colour  scales  are  used  to  illustrate  the  details  of  the  kernels. 
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Direct  P  arrival 
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Figure  3.  The  waveform  cross-correlation  measurements  of  the  traveltime  and  amplitude  anomalies  caused  by  the  cylindrical-shaped  wave  speed  perturbation 
are  compared  to  the  predicted  values  from  the  sensitivity  kernels.  The  horizontal  axes  are  traveltime  or  amplitude  anomaly,  and  the  vertical  axes  are  the  depth 
of  the  centre  of  the  cylindrical-shaped  wave  speed  perturbation  along  the  vertical  line  between  the  source  and  receiver  in  Fig.  1 .  The  left-hand  panels  are  for 
the  direct  P  arrival.  The  right-hand  panels  are  for  the  direct  S  arrival  when  a  double  couple  source  is  used.  The  top  row  is  for  the  model  with  only  P  wave 
speed  perturbations  (a  square  cosine  function  with  the  3  per  cent  maximum  at  the  centre  of  the  perturbation  and  zeros  at  the  boundaries),  the  middle  row 
for  the  model  with  only  S  wave  speed  perturbations  of  the  same  magnitude,  and  the  bottom  row  for  the  model  with  both  P  and  S  wave  speed  perturbations 
(3  per  cent  uniform). 


Normalized  Distance  (L//p) 

Figure  4.  The  relative  magnitudes  of  the  integrated  absolute  P  traveltime 
(asterisk)  and  amplitude  (triangle)  sensitivities  to  S  and  P  wave  speed  pertur¬ 
bations  vary  with  the  normalized  distance  between  the  source  and  receiver. 
The  P  and  S  wave  speeds  of  the  medium  and  the  source  time  function  are  the 
same  as  in  the  numerical  experiment  (Fig.  1).  At  one  wavelength  distance, 
the  ratios  are  110  per  cent  for  traveltime  and  540  per  cent  for  amplitude  (off 
scale). 


in  Dahlen  et  al.  2000).  This  causes  a  noticeably  wider  "doughnut’ 
hole  in  Kp  compared  to  that  of  K£  . 

Keeping  the  far-field  term  on  the  source  side  but  including  also 
the  near-  and  mid-field  terms  on  the  receiver  side  results  in  the  P 
traveltime  sensitivity  to  S  wave  speed  perturbations  along  the  ray 
path,  with  the  most  significant  change  on  the  receiver  side  (Fig.  6). 
The  change  of  Kpp  from  the  zero  sensitivity  along  the  ray  path  in 
the  case  with  only  the  far-field  term  (Fig.  5)  to  the  local  maximum 
sensitivity  in  Fig.  6  is  attributed  to  two  factors.  First,  the  scattering 
coefficient  for  a  mid-field  scattered  P  wave  due  to  an  S  wave  speed 
perturbation  is  the  maximum  along  the  ray  path  (Aki  &  Richards 
2002;  Favier  el  al.  2004).  Second,  the  mid-field  term  of  a  scattered 
wave  is  proportional  to  the  time  derivative  of  the  reference  displace¬ 
ment  at  the  receiver  (Aki  &  Richards  2002).  For  such  a  perturbed 
waveform,  the  delay  time  calculated  by  eq.  (Al)  in  Appendix  A 
is  the  maximum  when  there  is  no  offset  in  the  arrival  times  of  the 
reference  and  scattered  waves  (in  other  word,  when  the  scatter  is  on 
the  ray  path). 

Replacing  the  far-field  term  on  the  source  side  with  the  near-  and 
mid-field  terms  while  keeping  the  full  terms  on  the  receiver  side 
reveals  that  the  near-  and  mid-field  terms  on  the  source  side  also 
contribute  to  the  P  traveltime  sensitivity  to  S  wave  speed  perturba¬ 
tions  along  and  near  the  ray  path,  with  a  stronger  sensitivity  near  the 
source  (Fig.  7). 
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Figure  5.  Far-field  P  traveltime  and  amplitude  sensitivity  kernels  to  P  (top  panels)  and  S  (bottom  panels)  wave  speed  perturbations  in  the  vertical  plane 
containing  the  source  and  receiver.  The  units  of  the  sensitivities  are  the  same  as  in  Fig.  2. 


Figure  6.  The  P  traveltime  and  amplitude  sensitivity  kernels  to  P  (top  panels)  and  S'  (bottom  panels)  wave  speed  perturbations  are  calculated  with  the 
contributions  from  the  far-held  term  on  the  source  side  and  all  the  held  terms  on  the  receiver  side.  The  vertical  cross-section  contains  the  source  and  receiver 
and  the  units  of  the  sensitivities  are  the  same  as  in  Fig.  2. 


We  note  the  sensitivity  to  S  wave  speed  perturbations  oscillates 
from  negative  to  positive  values  away  from  the  ray  path  (Figs  2,  3 
and  6).  Unlike  the  sensitivity  to  P  wave  speed  perturbations,  inte¬ 
gration  of  the  sensitivity  to  S  wave  speed  perturbations  in  a  plane 
perpendicular  to  the  ray  path  causes  the  cancellation  of  the  positive 
and  negative  values.  A  uniform  S  wave  speed  perturbation  in  a  plane 
perpendicular  to  the  ray  path  yields  no  travel  time  anomaly  in  the  P 
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arrival.  This  means  that  S  wave  speed  perturbations  that  are  much 
larger  than  the  width  of  the  Fresnel  zone  (~  Isffd  when  /  >  1, 
where  /  is  the  distance  to  the  receiver  and  \  the  wavelength,  Dahlen 
et  al.  2000)  contribute  little  to  P  waveform  perturbations.  Integrat¬ 
ing  P  delay  times  caused  by  the  near-,  mid-  and  far-field  terms  in 
planes  perpendicular  to  the  ray  path,  Favier  et  al.  (2004)  noted  that 
the  far-field  term  becomes  predominant  at  distances  greater  than 
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Figure  7.  The  P  traveltime  and  amplitude  sensitivity  kernels  to  P  (top  panels)  and  S  (bottom  panels)  wave  speed  perturbations  are  calculated  with  the 
contributions  from  the  near-  and  mid-field  terms  on  the  source  side  and  all  the  field  terms  on  the  receiver  side.  The  vertical  cross-section  contains  the  source 
and  receiver.  The  units  of  the  sensitivities  are  the  same  as  in  Fig.  2. 


Figure  8.  A  cartoon  illustrating  wave  propagation  due  to  a  scatter.  The  star 
represents  a  source,  the  triangle  a  receiver,  and  the  circle  a  scatter.  If  only 
the  P  wave  speed  is  perturbed,  only  P  wave  is  scattered  by  the  scatter  (the 
upper  figure).  But  both  P  and  S  waves  are  scattered  if  the  S  wave  speed  is 
perturbed  (the  lower  figure). 

l/20th  wavelength  from  the  receiver  and  concluded  that  the  mid- 
field  term  cannot  be  neglected  at  distances  within  one  wavelength 
from  the  receiver.  However,  for  heterogeneities  on  the  scale  of,  and 
smaller  than,  the  width  of  the  sensitivity  kernels  (in  other  words  in 
places  where  the  finite-frequency  kernels  matter  for  tomography), 
the  absolute  or  maximum  sensitivities  are  more  meaningful  values 
to  assess  the  relative  importance  of  the  finite-frequency  kernels  to  P 
and  5  wave  speed  perturbations.  Fig.  4  shows  that  the  contributions 
of  5  wave  speed  perturbations  to  P  traveltime  and  amplitude  anoma¬ 


lies  are  substantial  at  distances  within  several  tens  of  wavelengths 
from  the  source  and  receiver.  So  this  cross-dependence  of  P  trav¬ 
eltime  and  amplitude  on  5  wave  speed  is  a  typical  finite-frequency 
effect,  important  where  the  scale  of  the  velocity  heterogeneity  is 
comparable  or  smaller  than  the  Fresnel  zone  of  the  wave. 

Our  results  differ  from  the  common  practice  in  tomographic  in¬ 
version,  in  which  P-wave  traveltimes  depend  only  on  P  wave  speed 
and  5-wave  traveltimes  only  on  5  wave  speed.  This  common  prac¬ 
tice  may  stem  from  two  well-known  theoretical  derivations:  the 
Helmholtz  poterntial  decomposition  in  a  homogenous  isotropic  me- 
diaum,  and  ray  theory,  in  which  the  spatial  gradients  of  the  elastic 
properties  are  neglected  under  the  high-frequency  approximation 
(e.g.  Dahlen  &  Tromp  1998;  Aki  &  Richards  2002).  The  effects  of 
5  wave  speed  perturbations  on  P  waveforms  presented  in  this  pa¬ 
per  result  from  the  terms  neglected  by  ray  theory.  For  general  elastic 
media,  the  spatial  gradients  of  the  elastic  properties  result  in  a  scalar 
(P)  wave  potential  that  explicitly  depends  on  the  spatial  gradient  of 
the  5  wave  speed  (Appendix  B).  The  finite-frequency  kernels  quan¬ 
titatively  reflect  that  there  is  only  scattered  P  wave  when  only  P  wave 
speed  is  perturbed.  But  an  5  wave  speed  perturbation  generates  both 
scattered  P  and  5  wave  (Fig.  8). 

Although  the  reference  model  used  in  this  paper  is  an  isotropic  and 
homogeneous  medium,  the  spatial  gradient  in  the  elastic  properties 
must  be  considered  once  a  wave  speed  perturbation  is  introduced 
into  the  medium.  The  numerical  approach  shown  in  Fig.  2  can  be 
applied  to  any  general  inhomogeneous  media,  in  which  the  finite- 
frequency  phenomena  documented  in  this  paper  should  also  exist. 


CONCLUSIONS 

For  finite-frequency  waves,  5  wave  speed  perturbations  may  have 
significant  effects  on  P  waveforms.  Whether  this  cross-dependence 
is  small  enough  to  be  neglected  depends  on  the  scale  of  wave  speed 
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heterogeneities  in  the  Earth  and  achievable  resolution  in  specific 
studies.  At  short  source-receiver  distances  (Fig.  4)  and  near  the 
source  and  receiver,  the  cross-dependence  of  P  waveforms  on  S 
wave  speed  is  substantial.  This  has  important  implications  for  seis¬ 
mic  tomography,  particularly  for  regional  teleseismic  tomography 
and  local  earthquake  tomography,  in  which  the  volume  of  interest 
is  often  within  several  10s  of  wavelength  from  receivers  or  both 
sources  and  receivers.  Even  for  regions  far  from  sources  and  re¬ 
ceivers,  neglecting  this  cross-dependence  means  that  the  estimate 
of  the  P  wave  speed  near  the  source  (e.g.  subduction  zone)  and  re¬ 
ceiver  is  likely  biased  and  such  errors  may  propagate  into  the  rest 
of  the  model  in  tomographic  inversions.  When  the  goal  of  tomo¬ 
graphic  inversion  is  to  place  constraints  on  thermal  and  composi¬ 
tional  variations  in  the  crust  and  mantle,  the  cross-dependence  of  P 
waveforms  on  S  wave  speed  variations  documented  above  becomes 
important  and  requires  a  joint  solution  of  both  the  P  and  S  wave  speed 
structures. 
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APPENDIX  A:  CALCULATION  OF  THE  SENSITIVITIES  BY  THE  SCATTERING 
INTEGRAL  METHOD 


We  define  the  traveltime  and  amplitude  anomalies,  8zp  and  8Aq ,  respectively,  as  (Tanimoto  1995;  Dahlen  et  al.  2000;  Zhao  et  al.  2005): 


P  ui(t)8u,(t)dt 
/,;  \Mt)\2dt 


(Al) 


ft ui(t)8ui(t)dt 

|S/(0l2d? 


(A2) 


where  U/(t)  is  the  displacement  calculated  from  the  reference  model  in  the  direction  e;,  a  unit  vector;  8u;(t)  —  ut(t)  —  Ui(t);  «/(f);  is  the 
displacement  for  the  perturbed  velocity  model;  a  dot  represent  the  time  derivative;  t  j  and  t2  are  the  upper  and  lower  limit  of  the  time  window  of 
the  arrival.  The  amplitude  anomaly  defined  here  is  different  from  that  in  Zhao  et  al.  (2005 )  by  a  constant.  Following  the  same  algebra  in  Zhao  et  al. 
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(2005),  we  have: 

1  f'2 

Ka,p=—  /  2paui(rR,t;rs ) 
n  Jt, 


1  f'2 

Ka,P=  —  2  pam(rR,t;rs) 

(2;  Jn 


If'2 

Kft,P=~^  2pfiui(rR,  f;  rs) 
n  Jh 


e,  •  [(V  •  Gr)(V  •  ii)]drdt 

e,  ■  [(V  •  Gr)(V  ■  ii)]drd? 

e,  ■  {((VG213)  :  [(VS)  +  (Vn)r]  -  2(V  •  Gr)(V  •  zi))drdf 


(A3) 

(A4) 


(A5) 


1  c{ 2  ~  c° 

=  “77  /  2pf3u,(rR,  t;rs )  / 
Jn  J-c 


e,  ■  {((VG213)  :  [(VS)  +  (VS)r]  -  2(V  •  Gr)(V  •  £)}drd/, 


(A6) 


where  the  G  is  the  Green’s  tensor.  The  symbol  ()213  represents  the  transposition  of  the  first  and  the  second  indices  of  a  third-order  tensor.  Pi 
and  Qi  are  the  normalization  factor  in  the  direction  e/,  which  are: 


(A7) 


(A8) 


For  the  direct  P  wave,  for  example,  the  sensitivity  kernels  to  P  and  S  wave  speed  perturbations,  Kpa  p  and  Kp  ,  are  calculated  from  eq.  (A3) 
and  (A5),  respectively,  with  a  time  window  (G,  t2 )  that  contains  the  P  arrival. 


APPENDIX  B:  WAVE  EQUATION  REPRESENTED  BY  THE  HELMHOLTZ  POTENTIAL 

The  cross-dependence  of  P  waveforms  to  5  wave  speed  perturbations  and  the  lack  of  effects  of  P  wave  speed  perturbations  to  S  waveforms 
can  be  understood  from  the  elastodynamic  equation, 

pu  =  V  ■  a  ,  (A9) 

where  p  is  density,  u  displacement,  and  a  stress.  For  simplicity,  the  body  force  term  is  omitted  in  the  equation.  The  displacement  can  be 
represented  by  the  Helmholtz  scalar  potential  of  P  wave,  cj>,  and  the  vector  potential  of  S  wave,  i/r  (e.g.  Aki  &  Richards  2002), 

u  =  V<p  +  V  x  i/f  with  V  ■  \[r  =  0.  (A10) 

Using  the  linear,  isotropic  stress-strain  relationship,  taking  the  divergence  and  curl  of  eq.  (A9),  and  assuming  a  constant  density,  we  have, 

=  V2(orV20)  +  2V  •  [VV  x  f  ■  VP2  +  VVtp  ■  V/S2  -  V£2V20]  (All) 

V  x  V  x  Vf  =  —V  x  V  x  [^2V  xVxf]  +  2Vx  [VV  x  f  ■  V^2  +  VV0  •  V,S2  -  V/S  2V20],  (A12) 

While  eq.  (A1 1)  shows  that  the  P- wave  scalar  potential  0  is  explicitly  related  to  the  P  wave  speed,  a ,  and  the  spatial  gradient  of  the  S  wave 
speed,  /S,  the  S-wave  vector  potential  i //  is  related  only  to  /S  in  eq.  (A12).  Although  eq.  (A12)  includes  two  terms  related  to  0,  they  have  no 
effect  on  \[r,  where  the  spatial  gradient  in  S  wave  speed  is  zero  or  when  \[r  and  cj>  are  separated  in  the  time  domain. 
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Shen,  Y.,  W.  Zhang,  Y.  Ren,  and  X.  Liang,  Imaging  the  crust  and  upper  mantle  in  the 
southern  and  southeastern  Tibetan  plateau:  A  three-dimensional  full-wavefield  approach, 
Eos  Trans.  AGU,  89(53),  Fall  Meet.  Suppl.,  abstract  T1  IE-06  Invited,  2008b. 

Finite-frequency,  teleseismic  traveltime  tomography  based  on  body-wave  ray  theory  in 
conjunction  with  the  Born  approximation  has  shown  wave-speed  anomalies  that  are 
indicative  of  the  delamination  of  the  mantle  lithosphere  and  possibly  part  of  the  lower 
crust  as  a  likely  cause  for  the  north-south-tending  rifts  in  the  southern  and  southeastern 
Tibetan  plateau.  The  images  also  show  evidence  that  contradicts  the  notion  of  substantial 
under-thrusting  of  the  Indian  mantle  lithosphere  beneath  southern  and  southeastern  Tibet. 
To  further  improve  resolution  and  verify  the  interpretations,  we  carry  out  a  three- 
dimensional  (3D)  full-wavefield  finite-frequency  tomography  study  of  the  area, 
incorporating  regional  earthquakes  and  ambient  signals  recorded  by  the  permanent  and 
portable  broadband  seismic  stations.  The  new  approach  accounts  for  complex  wave 
propagation,  enables  fuller  utilization  of  an  arrival  on  all  three  components  of  seismic 
records,  and  linearizes  the  inverse  problem  by  iteratively  updating  the  3D  reference 
model.  An  important  benefit  of  this  physically  realistic  modeling  of  full  wavefield  is  the 
consistency  of  the  system  of  equations  in  inversion,  which  is  particularly  important  for 
the  integration  of  different  types  of  observations  and  physical  properties.  We  calculate 
synthetic  waveforms  using  a  non-staggered  grid  finite-difference  code  in  a  polar- 
spherical  coordinate  system.  Frequency-dependent  phase  and  amplitude  anomalies  are 
measured  by  the  cross-correlation  between  synthetic  and  recorded  seismograms.  The 
finite-frequency  structural  sensitivity  kernels  of  these  measurements  are  calculated  by  the 
scattering-integral  method.  The  final  joint  3D  P-  and  S-speed  structural  model  is  obtained 
through  a  non-linear  iterative  process.  The  results  will  be  used  to  address  the  questions 
about  the  under-thrusting  of  the  Indian  mantle  lithosphere,  the  origin  of  the  north-south- 
trending  rifts,  crustal  and  mantle  lithosphere  delamination,  and  their  implications  for  the 
rise  of  the  Tibetan  plateau. 
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We  carry  out  a  full-wave  finite-frequency  tomography  of  the  Tibetan  plateau  using  the 
scattering-integral  method.  We  use  about  110  stations  from  the  PASSCAF  experiments  in 
the  region,  which  include  the  HIMNT,  MIT-China  and  Namche-Barwa  experiments.  We 
start  by  re-locating  the  earthquakes  and  re-detennining  the  focal  mechanisms  using  3-D 
strain  Green  Tensors.  The  Harvard  CMT  solutions,  if  available,  are  used  as  initial 
solutions  of  the  iteration.  We  then  calculate  the  synthetic  wave  field  using  the  updated 
source  locations  and  focal  mechanisms  in  a  3D  reference  model.  Finite-frequency 
traveltimes  measured  by  cross-correlating  observed  and  synthetic  waveforms  are  inverted 
for  wave-speed  perturbations  to  the  3D  reference  model.  We  iteratively  redo  the  two 
steps  to  improve  the  quality  of  the  tomographic  images.  After  two  iterations,  the  model 
reveals  strong  velocity  heterogeneities,  with  local  changes  of  ±10%  relative  to  the  initial 
reference  model  (CUB  2.0).  It  shows  low  velocity  anomalies  between  the  Bangong- 
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Nujiang  Suture  and  the  Indus-Tsangpo  Suture  for  both  Vp  and  Vs  velocities  in  a  depth 
range  of  30km  to  100km.  This  result  implies  that  the  Indian  mantle  lithosphere  does  not 
underthrust  beyond  the  Indus-Tsangpo  Suture  in  eastern  Tibet. 
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